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INTRODUCTION: 


This  final  technical  report  comprises  five  sections  representing  different  facets  of 
work  undertaken  at  the  Research  School  of  Earth  Sciences,  Australian  National 
University  on  the  problems  of  seismic  wave  propagation  and  inference  about  three- 
dimensional  structure  from  regional  to  teleseismic  scales. 

Part  I:  Upper  mantle  beneath  Australia  from  {portable  array  deployments. 

R.  van  der  Hilst  &  B.L.N.  Kennett 

This  contribution  brings  together  the  range  of  studies  which  have  been  made  on  the 
crustal  and  upper  mantle  structure  beneath  Australia  exploiting  portable  station 
deployments,  notably  the  Skippy  array.  A  detailed  three-dimensional  model  of  shear 
wave  speed  has  now  been  produced  for  two-thirds  of  the  continent  and  the  oceanic 
zone  to  the  east,  in  less  than  four  years  from  the  commencement  of  field  recording. 

The  integration  of  receiver  based  studies  including  P  wave  delay-time  tomography 
and  crustal  receiver  functions  with  the  model  derived  from  surface  wave  inversion 
provides  an  excellent  characterisation  of  a  continent.  Such  a  representation  can  be 
used  to  calibrate  both  regional  and  teleseismic  information. 

Part  II:  Upper  mantle  heterogeneities  in  the  Indian  Ocean  from  waveform  inversion. 
E.  Debayle  &  }.}.  Levique 

This  section  represents  work  begun  at  the  University  of  Strasbourg  in  France  and 
completed  at  the  Australian  National  University.  The  application  of  joint  inversion 
of  Rayleigh  and  Love  wave  trains  to  extract  information  on  both  velocity 
heterogeneity  and  the  pattern  of  upper  mantle  anisotropy  is  of  particular  interest  in  a 
continental  setting.  A  comparable  inversion  is  underway  for  the  SKIPPY  data  across 
the  Australian  continent. 

Part  III:  Guided  waves  in  3-dimensional  structures 
B.L.N.  Kennett 

Current  treatments  of  surface  wave  propagation  in  three-dimensionally  varying 
media  depend  on  perturbation  techniques.  For  2-D  models  it  is  possible  to  make  a 
full  development  which  describes  the  evolution  of  the  wavefield  in  terms  of  coupling 
between  the  surface  wave  modes  of  a  reference  model.  The  extension  to  3-D  models 
requires  simultaneous  consideration  of  mode-coupling  and  transfer  between  plane 
wave  components  so  that  phenomena  such  as  surface  wave  refraction  can  be 
described.  The  theory  is  quite  complex  but  for  periods  in  the  range  25-100  seconds 
should  provide  a  suitable  basis  fora  computational  development  using  a  limited 
number  of  modes. 


Part  IV:  A  regionalized  upper-mantle  (RUM)  model. 

O.  Gudmundsson  &  M.  Sambridge 

Models  of  three-dimensional  structure  in  the  upper-mantle  derived  from 
conventional  travel-time  tomography  have  coverage  restricted  to  the  vicinity  of 
seismic  sources  and  receivers.  However,  we  can  extend  the  information  by 
representing  the  major  features  around  the  globe  using  a  relatively  detailed  tectonic 
regionalization  and  then  undertaking  an  inversion  for  the  best  fitting  upper  mantle 
for  each  such  tectonic  class.  The  resulting  model  provides  a  means  of  estimating  the 
properties  of  the  upper  mantle  e.g.  teleseismic  residuals  even  in  zones  which  have 
not  yet  been  characterised.  Such  a  procedure  can  therefore  provide  a  first  step 
towards  the  calibration  of  paths  to  compensate  for  three-dimensional  structure  in  the 
process  of  event  location. 
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SUMMARY 

The  distribution  of  earthquakes  at  regional  distances  around  Australia  is 
particularly  favourable  for  these  natural  events  to  be  used  as  probes  into  the 
seismic  structure  of  the  lithosphere  and  upper  mantle.  The  distribution  of 
permanent  seismic  stations  is  too  sparse  for  detailed  continent-wide  seismic 
imaging  and  most  information  used  in  our  studies  has  come  from  deployments 
of  portable  seismic  recorders. 

Such  portable  instruments  were  initially  used  in  northern  Australia  to 
define  the  radial  variations  in  P  and  S  velocity  and  attenuation  using  a 
combination  of  short-period  and  broadband  observations,  and  to  investigate 
seismic  anisotropy  from  shear-wave  splitting  observed  in  refracted  S  waves 
propagating  through  the  upper  mantle. 

The  whole  Australian  continent  has  now  been  covered  in  the  Skippy  experiment 
from  1993-1996  in  which  a  sequence  of  deployments  of  up  to  12  recorders  at 
a  time  have  been  used  to  synthesise  a  continental  scale  array  of  broad-band 
instruments.  The  major  objectives  of  this  project  are  to  delineate  the 
three-dimensional  structure  of  the  Australian  lithosphere  and  underlying 
mantle  by  tomographic  inversions  for  lateral  variations  in  seismic  wavespeed 
and  site  specific  studies  to  map  crustal  thickness  (receiver  functions)  and 
seismic  anisotropy  (shear  wave  splitting). 

Surface  wave  studies  have  begun  to  reveal  the  three-dimensional  variations  in 
shear  wave  structure  beneath  the  continent  by  exploiting  the  records  from 
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portable  broadband  stations.  Dense  data  coverage  enables  the  imaging  of 
wavespeed  variations  on  length  scales  of  250  km  and  larger  beneath  most 
of  central  and  eastern  Australia. 

Results  of  the  waveform  inversion  suggest  that  the  eastern  edge  of  the 
Proterozoic  shields  of  central  Australia  is  a  complex  three-dimensional  surface 
which  does  not  have  a  simple  relation  to  the  conventional  Tasman  Line 
marking  the  separation  of  Precambrian  and  Phanerozoic  outcrop.  In  the 
shallow  lithosphere  (80  km  depth)  the  most  pronounced  lateral  contrast  in 
seismic  wave  speed  occurs  significantly  further  east  than  the  conventional 
Tasman  Line.  However,  in  northern  Australia  the  eastern  boundary  of  exposed 
Precambrian  basement  coincides  with  a  transition  from  moderately  high 
wavespeeds  to  fast  wave  propagation  that  is  particularly  prominent  at  a  depth 
of  about  140  km.  At  larger  depth  still  (200  km)  the  wavespeed  gradient  seems 
to  occur  significantly  further  west  (near  the  western  margin  of  the  Eromanga 
basin).  Beneath  most  of  eastern  Australia  the  fast  seismic  ‘lid’  does  not  extend 
to  depths  larger  than  about  100  km,  and  there  is  a  pronounced  low  velocity 
zone  between  about  100  and  200  km  depth.  The  surface  wave  data  do  not 
indicate  a  slow  wavespeed  channel  beneath  central  Australia.  Instead,  the  fast 
lid  extends  to  a  depth  of  at  least  250  km,  and  locally  to  depths  in  excess  of  350 
km.  The  north  Australian  craton  is  seismically  well  defined  but  does  not  seem 
to  continue  into  the  Kimberley  block.  In  the  shallow  lithosphere,  the  region 
influenced  by  the  late  Palaeozoic  Alice  Springs  orogeny  (Amadeus  basin  and 
Musgrave  block)  is  marked  by  significantly  slower  seismic  wave  propagation 
than  the  Proterozoic  cratons  to  the  north  and  south. 

Analysis  of  data  from  individual  stations  has  produced  a  set  of  shear  velocity 
profiles  for  the  crust  and  uppermost  mantle  beneath  Skippy  stations  which 
provide  a  useful  complement  to  the  information  at  depth  from  the  surface 
wave  studies.  The  crust-mantle  boundary  is  deep  (38-44  km)  and  mostly 
transitional  in  character  along  the  axis  of  the  eastern  fold  belt,  but  is  relatively 
sharp  and  shallower  (30-36  km)  near  the  boundary  between  Precambrian  and 
Phanerozoic  outcrop.  Shear  wave  splitting  results  from  the  portable  stations 
are  beginning  to  reveal  a  complex  pattern  of  seismic  anisotropy  that  requires 
further  analysis  of  both  body  and  surface  wave  data. 
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1  INTRODUCTION 

The  surficial  geology  of  the  Australian  continent  is  composed  of  an  assemblage  of 
crustal  blocks  that  can  be  broadly  grouped  into  the  Precambrian  western  and  central 
cratons  and  the  Phanerozoic  eastern  province  (figure  lA).  Structural  differences  in 
the  mantle  beneath  the  Precambrian  shield  and  eastern  Australia  have  previously 
been  inferred  from  limited  observations  of  surface  wave  dispersion  (cf.  Muirhead  & 
Drummond,  1991;  Denham,  1991),  and  teleseismic  travel-time  residuals  (Drummond 
et  al.,  1989).  These  results  have  been  confirmed  and  extended  in  recent  studies  using 
deployments  of  portable  broad-band  instruments  across  the  whole  continent  (the  Skippy 
experiment.  Van  der  Hilst  et  al.,  1994),  and  it  is  now  possible  to  derive  three-dimensional 
models  of  the  P  and  S  wavespeeds  beneath  the  Australasian  region  to  depths  of  400  km 
or  more. 

The  extensive  earthquake  activity  in  the  seismic  belt  that  runs  through  Indonesia, 
New  Guinea  and  its  offshore  islands,  Vanuam,  Fiji  and  the  Tonga-Kermadec  zone 
provides  a  wide  range  of  natural  sources  which  can  be  used  to  constrain  seismic 
structure  in  the  lithosphere,  asthenosphere,  and  the  transition  zone  beneath.  There  is 
a  rather  limited  number  of  high-quality,  permanent  seismological  stations  in  Australia 
so  that  for  high  resolution  studies  of  seismic  strucmre  these  observatories  need  to  be 
supplemented  by  the  installation  of  portable  recorders  (figure  IB). 

Building  on  the  experience  from  array  deployments  in  northern  Australia  the 
Research  School  of  Earth  Sciences  (RSES)  started  in  1993  a  nationwide  observational 
seismology  project,  the  Skippy  experiment,  which  uses  a  sequence  of  array  deployments 
to  synthesise  a  continental  scale  array  of  broad-band  instruments.  Major  objectives  of 
this  project  are  to  exploit  different  classes  of  seismological  data  for  the  delineation  of 
the  three-dimensional  structure  of  the  Australian  lithosphere  and  underlying  mantle 
by  tomographic  inversions  for  lateral  variations  in  seismic  wavespeed  and  site  specific 
studies  to  extract  crustal  Information  through  the  construction  of  receiver  functions, 
and  seismic  anisotropy  from  analysis  of  shear  wave  splittinc. 

For  instruments  in  northern  Austraha,  refracted  arrivals  for  both  P  and  S  waves  from 
the  upper  mantle  can  be  used  to  determine  upper  mantle  structure.  At  greater  distances 
from  the  sources  multiple  S  arrivals  and  surface  waves  provide  a  powerful  probe  for 
three-dimensional  structure  using  information  on  many  crossing  propagation  paths. 
For  this  class  of  arrival  the  energy  is  largely  directed  horizontally  and  the  waveforms 
are  most  sensitive  to  S  wave  velocities.  The  information  from  regional  earthquakes  can 
be  supplemented  with  the  arrival  times  for  teleseismic  waves  in  delay-time  tomography 
for  P  arrivals  where  the  paths  through  the  mantle  are  relatively  steep.  In  addition, 
information  from  distant  events  can  be  used  to  constrain  crustal  and  upper  mantle 
structure  through  the  analysis  of  converted  waves  and  reverberations. 

The  combination  of  many  different  classes  of  seismological  information  has  begun  to 
reveal  the  complex  nature  of  lithospheric  and  mantle  structure  and  to  shed  light  on  the 
way  in  which  the  Australian  continent  may  have  been  assembled.  In  this  paper  we  briefly 
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review  results  of  the  array  experiments  in  northern  Australia,  from  both  short-period 
and  broadband  instruments,  and  present  the  latest  results  from  the  Skippy  project  for 
the  continent  as  a  whole. 


2  P  AND  S  VELOCITY  PROFILES  IN  THE  UPPER  MANTLE 

The  events  to  the  north  of  Australia  provide  a  convenient  set  of  energy  sources  for 
studies  of  the  upper  mantle  and  it  has  proved  possible  to  constrain  the  major  features 
of  the  mantle  velocity  profile  as  well  as  the  attenuation  distribution  with  depth  for  both 
P  and  S  waves. 

2.1  Short-period  studies 

Much  of  the  information  on  mantle  structure  has  come  from  deployments  of  portable 
instruments  with  short-period  seismometers.  Hales,  Muirhead  &  Rynn  (1980)  carried  out 
a  travel  time  analysis  for  Indonesian  earthquakes  that  occurred  at  a  variety  of  depths 
and  were  recorded  at  a  number  of  portable  stations  in  northern  Australia.  The  resulting 
model  is  rather  complex  with  many  small  discontinuities  and  low  velocity  zones,  which 
may  reflect  the  mapping  of  three-dimensional  structure  into  a  one-dimensional  profile.  A 
subsequent  re-interpretation  by  Leven  (1985),  using  comparisons  between  observed  and 
synthetic  seismograms,  leads  to  somewhat  simplified  structure  but  retains  a  prominent 
velocity  contrast  near  210  km  depth. 

In  the  period  1985-1987,  RSES  carried  out  a  sequence  of  experiments  using 
short-period  vertical  seismometers  in  northern  Australia  (figure  IB)  to  record  the  natural 
seismicity  in  the  Indonesia/New  Guinea  region.  Many  of  the  results  were  summarised 
by  Dey  et  al  (1993)  who  present  composite  record  sections  of  upper  mantle  arrivals 
that  show  significant  variation  in  P  wave  velocity  structure  between  paths  for  events 
along  the  Flores  arc,  studied  by  Bowman  &  Kennett  (1991),  and  paths  to  events  in  New 
Guinea.  The  shallow  structure  has  to  be  inferred,  but  the  P  velocity  structure  is  well 
constrained  from  above  the  base  of  the  lithosphere  near  210  km  down  to  below  the 
410  km  discontinuity.  The  interpretation  confirms  the  need  for  a  P  velocity  contrast 
near  210  km  depth.  For  S  waves  the  corresponding  record  sections  show  a  clear  arrival 
associated  with  the  lithosphere  which  cannot  easily  be  traced  beyond  2000  km  but  no 
branches  associated  with  greater  depth. 

2.2  Broad-band  studies 

Modern  broadband  seismometers  provide  a  faithful  rendition  of  ground  motion  over 
a  wide  range  of  frequency  and  allow  the  full  exploitation  of  both  P  and  S  body  waves. 
A  broadband  sensor  has  been  operated  by  RSES  at  the  Warramunga  array  in  Northern 
Australia  since  late  in  1988  (station  code  WRA);  this  facility,  was  upgraded  in  1994  by 
IRIS  at  a  nearby  site  (WRAB).  Over  a  period  of  years  it  has  been  possible  to  build  up 
record  sections  covering  the  range  of  interest  for  the  upper  mantle  by  using  events  in 
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the  Indonesia/New  Guinea  earthquake  belt.  The  records  from  the  permanent  station 
have  been  augmented  by  portable  broadband  stations  deployed  at  distances  of  up  to 
300  km  from  WRA.  The  surface  conditions  in  this  region  are  such  that  good  results 
can  be  obtained  for  SV  waves  on  radial  component  records  after  rotation  to  the  great 
circle  path:  the  high  surface  velocities  lead  to  little  contamination  by  converted  P  waves. 
This  represents  a  considerable  benefit  over  previous  S  wave  studies  of  the  upper  mantle 
which  have  been  restricted  to  SH  waves  (Gudmundsson,  Kennett  &  Goody,  1994). 

The  Skippy  experiment  (Van  der  Hilst  et  al.,  1994)  has  emplaced  a  number  of  portable 
broadband  instruments  in  northern  Australia  at  a  suitable  distance  range  to  complement 
the  observations  at  the  Warramunga  array.  Figure  2  shows  a  composite  record  section 
covering  the  P  and  S  wave  components  returned  from  the  upper  mantle  for  events  in 
New  Guinea  recorded  at  Skippy  stations  in  Queensland.  This  section  has  been  constructed 
from  unflltered  vertical  component  records  from  five  shallow  events  and  clearly  displays 
the  benefit  of  broadband  recording.  The  onset  of  the  S  waves  shows  high  frequency 
behaviour  (greater  than  1  Hz)  out  to  2000  km,  but  beyond  this  distance  the  S  wave 
arrivals  have  a  significantly  lower  frequency  (0.2  Hz  at  3000  km)  and  this  is  also  seen 
for  later  arrivals  at  shorter  distance.  For  P  waves  the  loss  of  higher  frequencies  is  less 
pronounced.  An  enlargement  of  the  S  wave  arrivals  for  both  the  radial  and  tangential 
components  is  displayed  in  figure  3  using  unfiltered  broadband  records. 

The  change  in  frequency  content  for  S  waves  returned  from  greater  depth  has 
been  analysed  by  Gudmundsson,  Kennett  &  Goody  (1994)  to  determine  the  attenuation 
structure  with  depth  tmder  northern  Australia.  The  slope  of  the  spectral  ratio  between 
P  and  S  wave  arrivals  on  the  same  record  has  been  used  to  determine  the  differential 
attenuation  between  P  and  S.  This  differential  Information  can  be  interpreted  with 
a  knowledge  of  the  velocity  structure  and  requires  strong  attenuation  of  S  waves  in 
the  asthenosphere  between  210  km  and  the  410  km.  In  a  parallel  analysis  Kennett, 
Gudmundsson  &  Tong  (1994)  have  used  the  composite  record  sections,  together  with 
the  earlier  information  from  the  short  period  studies,  to  buUd  velocity  profiles  for  P  and 
S.  These  velocity  models  have  been  refined  by  comparison  of  observed  and  synthetic 
seismograms  including  the  influence  of  attenuation  (figure  4).  An  advantage  of  this 
study  is  that  both  P  and  S  velocity  profiles  are  determined  for  the  same  events  and  so 
the  P/S  velocity  ratio  can  be  well  determined,  which  is  particularly  useful  for  studies 
of  mantle  composition.  The  depth  variation  of  the  P/S  velocity  ratio  is  in  good  general 
agreement  with  the  results  for  the  shield  areas  of  north  America  obtained  by  combining 
the  P  velocity  profile  of  LeFevre  &  Helmberger  (1989)  with  the  S  wave  structure  of  Grand 
&  Helmberger  (1984). 

Comparison  of  the  radial  and  tangential  components  of  the  S  wavefield  at  WRA 
indicates  that  for  the  arrivals  from  the  upper  mantle  discontinuities  there  is  a 
systematically  earher  arrival  on  the  SH  component  by  more  than  a  second.  These 
indications  of  seismic  anisotropy  for  the  refracted  S  wave  arrivals  have  been  analysed 
by  Tong,  Gudmundsson  &  Kennett  (1994)  who  have  determined  the  direction  of  fast 
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propagation  and  time  shift  between  the  S  components  by  correlation  analysis.  The 
nature  of  the  observed  anisotropy  is  such  that  it  cannot  be  explained  by  structure  local 
to  the  receivers.  A  level  of  anisotropy  of  the  order  of  1  percent  in  both  the  lithosphere 
and  the  asthenosphere  beneath  would  explain  the  data  quite  well. 


3  THREE-DIMENSIONAL  STRUCTURE 

The  favourable  position  of  the  Australian  continent  relative  to  world  seismicity  can  be 
exploited  in  a  number  of  ways  to  obtain  information  on  the  three-dimensional  seismic 
structure  in  the  mantle.  This  can  be  done  by  means  of  tomographic  techniques  that  rely 
on  the  number  of  high-quality  crossing  paths  produced  by  many  source-receiver  pairs 
or  by  site  specific  studies  that  aim  to  define  the  subsurface  structure  beneath  individual 
stations  by  minimizing  the  influence  of  structure  elsewhere.  These  applications  provide 
complimentary  information  about  Earth’s  structure,  which  are  being  merged  to  obtain 
a  more  complete  description  of  the  Australian  mantle. 

The  configuration  of  the  seismicity  around  the  Australian  continent  is  very  favourable 
for  tomographic  methods  since  with  a  suitable  distribution  of  receivers  a  dense  and  even 
data  coverage  of  the  continent  can  be  achieved.  RSES  has  just  completed  a  three  and 
a  half  year  field  program  to  install  some  60  portable  broadband  stations  across  the 
entire  Australian  continent  (see  figure  IB).  This  project  uses  a  set  of  up  to  12  portable 
instruments  which  occupy  sites  for  at  least  5  months  at  a  time  before  being  moved  to 
a  new  location,  so  that  a  full  continental  array  can  be  synthesised  (van  der  HUst  et  al 
1994,  Kennett  &  van  der  Hilst  1996).  The  mobility  of  the  arrays  have  led  to  the  name 
Skippy  (nickname  for  kangaroo)  for  the  whole  project.  The  deployments  commenced 
in  May  1993  with  8  stations  in  Queensland  and  the  coverage  of  the  whole  continent 
was  completed  at  the  end  of  September  1996.  The  data  set  from  the  Skippy  stations  is 
supplemented  by  records  of  suitable  events  from  the  permanent  broadband  stations.  In 
addition  to  providing  data  from  sites  not  occupied  by  Skippy  stations,  the  records  from 
the  observatory  instruments  are  important  in  constraining  upper  mantle  structure  from 
fundamental  mode  surface  wave  since  at  frequencies  less  than  10  mHz  they  are  often 
of  better  quality  than  the  records  of  the  portable  stations. 


3.1  Surface  wave  studies 

The  records  from  the  broadband  stations  have  been  used  in  a  number  of  different 
studies.  A  major  objective  is  the  delineation  of  lithospheric  and  mantle  structure  using 
waveform  tomography  for  the  shear  wave  and  surface  wave  portion  of  the  seismogram. 
So  far  the  analysis  has  been  based  on  the  partitioned  waveform  inversion  technique 
introduced  by  Nolet  (1990).  A  nonlinear  optimisation  is  used  to  find  a  stratified  model 
which  gives  the  best  fit  to  an  observed  seismogram,  which  should  represent  the  average 
structure  along  the  great  circle  between  source  and  receiver.  The  assemblage  of  path 
averages  are  then  used  in  a  linear  inversion  to  recover  the  three-dimensional  shear  wave 
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structure.  The  tomographic  inversion  of  the  linear  constraints  on  the  three-dimensional 
velocity  structure  uses  a  block  structure  of  approximately  equal  area  in  latimde  and 
longitude  coupled  into  a  a  set  of  mostly  triangular  basis  functions  in  the  vertical 
direction.  Both  model  norm  and  gradient  damping  are  used  to  achieve  a  balance  between 
data  fit  and  smoothness  of  the  model.  The  details  of  the  partitioned  waveform  inversion 
method  are  given  by  Nolet  (1990),  and  its  application  to  Australian  data  is  presented  by 
Zielhuis  &  Van  der  HUst  (1996). 

The  non-linear  waveform  technique  assumes  the  independent  propagation  of 
individual  modes  along  a  great  circle  path.  The  neglect  of  inter-mode  coupling  may 
prevent  the  reliable  imaging  of  deep  mantle  structure  and  we  therefore  restrict  the 
present  discussion  to  the  upper  mantle.  Since  we  do  not  account  for  out-of-plane 
scattering  we  do  not  make  use  of  high-frequency  fundamental  mode  surface  waves  that 
are  sensitive  to  large  wavespeed  variations  in  the  shallow  lithosphere  (such  as  variations 
in  crustal  thickness). 

The  effective  use  of  the  non-linear  inversion  scheme  for  the  individual  paths  requires 
synthetic  seismograms  that  are  close  to  the  observed  waveforms,  and  in  turn  this 
requires  a  reference  model  that  is  close  to  the  averaged  ID  velocity  model  for  the 
wave  path.  Care  has  to  be  taken  to  account  for  the  best  averaged  crustal  model  for 
the  path  in  a  region  which  covers  both  continental  and  oceanic  lithosphere.  A  partial 
compensation  for  crustal  thickness  variations  is  made  by  using  different  mode  fdes  for 
source,  propagation  path,  and  receiver  structure  (Zielhuis  &  Nolet,  1994b,  Kennett  1995). 
The  influence  of  heterogeneity  is  rather  different  for  the  fundamental  and  higher  modes. 
For  instance,  at  frequencies  higher  than  about  30  mHz  the  fundamental  modes  are 
sensitive  to  variations  in  crustal  structure  in  general  and  the  ocean-continent  transition 
in  particular.  The  higher  modes  are  more  sensitive  to  structure  at  larger  depth  and  are 
less  influenced  by  scattering  due  to  strong  heterogeneity  in  the  shallow  lithosphere  and 
can  thus  be  interpreted  to  higher  frequency.  The  waveform  matching  was,  therefore, 
restricted  to  frequencies  of  up  to  25  mHz  for  the  fimdamental  mode  and  50  mHz  for  a 
time  window  comprising  covering  the  range  of  phase  velocities  expected  for  the  higher 
modes.  The  use  of  higher  modes,  which  are  well  excited  by  the  deep  earthquakes  in  the 
region  (Tonga,  New  Hebrides,  Java),  improves  resolution  of  structure  in  depth.  However, 
the  influence  of  the  heterogeneity  in  structure  means  that  it  is  not  always  possible  to 
fit  both  the  fundamental  mode  and  higher  mode  windows  with  the  same  structure.  In 
future  studies  this  linearization  problem  wiU  be  remedied  by  using  wavespeed  profiles 
and  associated  mode  files  of  reference  models  that  are  constructed  specifically  for  the 
region. 

Zielhuis  &  Van  der  Hilst  (1996)  discussed  results  of  the  application  of  this  partitioned 
waveform  inversion  scheme  to  the  data  from  stations  in  eastern  Australia  (arrays  SKI, 
SK2,  and  BAS,  Figure  IB).  These  data  are  consistent  with  a  pronounced  change  from  slow 
shear  wave  propagation  beneath  the  Coral  and  Tasman  Seas  and  the  coastal  regions 
of  continental  Australia  to  very  fast  wave  propagation  beneath  central  Australia.  The 
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thickness  of  the  high  wavespeed  lid  increases  from  about  80-100  km  in  easternmost 
Australia  to  at  least  250  km  beneath  the  central  shields.  The  transition  is  rather  sharp 
in  north  Queensland  but  is  manifest  as  a  multiple  boundary  further  south.  Based  on 
these  data  Zielhuis  &  Van  der  Hilst  argued  that  the  the  conventional  Tasman  line  (figure 
IB)  marking  the  limits  of  Precambrian  outcrop,  may  roughly  coincide  wth  a  contrast  in 
seismic  wave  speed  in  the  north,  but  that  high  wavespeed  lithosphere  extends  east  of 
the  Precambrian  basement  further  to  the  south.  However,  the  data  used  in  this  earlier 
smdy  does  not  constrain  the  southermnost  part  of  the  region  around  the  Tasman  Line, 
nor  does  it  reliably  image  possible  lateral  variations  in  wavespeed  in  the  upper  mantle 
beneath  central  Australia. 

The  model  presented  here  is  based  on  waveform  data  from  the  permanent 
observatories  in  the  region  and  from  45  Skippy  stations  covering  the  eastern  and 
central  parts  of  the  Australian  continent  (arrays  SK1,SK2,SK3,SK4  and  BAS  in  Figure  IB). 
Low-noise  records  were  selected  for  regional  earthquakes  for  which  a  Centroid  Moment 
Tensor  (CMT)  is  available  from  Harvard  University.  The  total  number  of  seismograms 
used  in  the  construction  of  the  three-dimensional  shear  velocity  mode  is  approximately 
1350  from  about  360  earthquakes.  Good  wave  path  coverage  is  achieved  for  eastern 
and  central  Australia  as  shown  in  figure  5  which  enables  the  extension  of  the  results  for 
eastern  Australia  by  Zielhuis  &  Van  der  Hilst  (1996). 

Results  of  the  inversion  for  three-dimensional  variations  in  shear  wavespeed  are 
presented  in  figure  6  for  depths  of  80,  140,  200  and  300  km.  On  wavelengths  of  1000 
km  and  larger  the  observations  are  in  good  agreement  with  inferences  from  global 
inversion  of  Rayleigh  wave  phase  velocities  (e.g.,  Trampert  &  Woodhouse  1995),  but  the 
data  coverage  provided  by  the  Skippy  arrays  allows  the  study  of  continental  strucmre  in 
much  more  detail.  The  resolution  length  is  about  300  km  for  most  of  the  area  discussed 
here. 

A  striking  feature  in  the  images  of  mantle  strucmre  to  about  200  km  depth  is  the  large 
difference  in  shear  wavespeed  between  the  part  of  Australia  east  of  140°E  and  the  rest 
of  the  continent  (Figures  6A,B).  The  wavespeed  gradients  are  predominantly  in  west-east 
direction  and  exhibit  peak-to-peak  arnplimde  variations  of  up  to  10  per  cent  relative  to 
the  reference  model  used  for  the  display.  The  patterns  generally  match  geological  trends 
observed  at  the  surface.  Beneath  easternmost  Australia  and  the  adjacent  oceanic  regions 
the  shear  wave  speed  is  relatively  low,  whereas  fast  wave  propagation  marks  the  central 
and  western  part  of  the  continent.  The  boundary  between  these  regions  is  well  defined 
at  80  and  140  km  depth.  The  transition  from  low  to  high  velocities  seems  to  be  rather 
sharp  in  northern  Queensland  (near  15°S,  142°E)  and  lies  close  to  the  Tasman  Line,  the 
presumed  eastern  edge  of  the  Precambrian  shields  based  on  geological  outcrop,  gravity 
and  magnetic  data  (Figure  IB).  Further  to  the  south  the  transition  from  slower  to  faster 
wavespeeds  appears  to  be  somewhat  more  complex  and  may  occur  over  at  least  two 
zones  of  rapid  Increase  in  wavespeed.  The  southern  part  of  the  Eromanga  basin  and 
the  region  associated  with  the  Lachlan  fold  belt  is  characterized  by  wavespeeds  that  are 
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intermediate  to  that  in  the  coastal  regions  and  the  central  shields  (figure  6B).  The  nature 
of  these  wavespeed  transitions  and,  in  particular,  the  complex  wavespeed  structure  near 
the  Broken  Hill  block  {-32°S,  140°E)  will  be  discussed  in  more  detail  in  a  separate  paper 
(Kennett  &  Van  der  Hilst,  in  preparation).  We  remark  that  south  of  130°S,  the  velocity 
contrast  associated  with  the  boundary  between  Phanerozoic  and  Proterozoic  basement 
is  significantly  better  constrained  than  in  the  study  by  Zielhuis  &  Van  der  HUst  (1996). 

Along  the  eastern  seaboard  of  the  Australian  continent  there  is  a  pronounced  low 
velocity  zone  between  100  and  200  km  in  depth  (figure  6);  the  minimum  wavespeed 
seems  to  occur  at  a  depth  of  140  km.  This  low  velocity  zone  was  recognised  in  earlier 
surface  wave  dispersion  studies  (see  Muirhead  &  Drummond  1991  for  a  review)  but  a 
striking  feature  in  our  3-D  models  is  the  level  of  lateral  variation  in  the  character  of  the 
zone.  The  low  velocity  zone  is  interrupted  near  30°S,  150°E  beneath  the  New  England 
fold  belt  and  the  higher  velocities  in  this  region  seem  to  be  almost  continuous  from  the 
surface  to  the  transition  zone.  The  main  zones  of  lowered  velocities  correlate  quite  well 
with  recent  vulcanlsm  and  high  heat  flow.  The  low  speed  anomalies  appear  to  extend 
to  greater  depth  beneath  the  Coral  Sea  and  beneath  southeastern  corner  of  Australia. 
Zielhuis  &  Van  der  Hilst  (1996)  discuss  the  upper  mantle  strucmre  beneath  easternmost 
Austraha  in  more  detail. 

In  the  upper  200  km  beneath  central  Australia  shear  wavespeed  is  high,  but  the 
images  suggest  significant  lateral  variations  that  begin  to  delineate  the  major  crustal 
elements.  The  zone  of  high  wavespeed  in  the  central  north  (Figures  6A-C)  coincides  with 
the  lateral  extent  of  the  North  Australian  craton,  and  the  high  wavespeeds  in  the  central 
south  coincide  with  the  Gawler  craton  and  the  Eucla  block  (Figures  6A,B).  At  80  km 
depth,  the  region  associated  with  the  Alice  Springs  orogeny  -  in  particular  the  Amadeus 
basin  and  the  Musgrave  block  -  is  characterized  by  slower  shear  wave  propagation  than 
in  the  adjacent  shields  (Figure  6A).  The  images  suggest  that  the  north  Australia  craton 
extends  northward  into  western  Papua  New  Guinea  and  eastern  Indonesia  (from  Timor 
to  Man  Jaya).  Interestingly,  the  Kimberley  block,  near  15°S,  125°Eis  delineated  by  slower 
wave  propagation  than  the  shield  further  to  the  east,  which  suggests  that  the  Kimberley 
basin  is  not  underlain  by  the  same  basement  as  the  Proterozoic  North  Australia  craton. 
The  change  to  slower  wave  propagation  occurs  across  the  Ord  basin  and  the  Halls  Creek 
shear  zone.  The  high  wavespeeds  mapped  in  the  south  west  may  well  be  related  to  the 
Archaean  cratons  but  they  are  not  yet  well  localized  owing  to  insufficient  path  coverage. 

The  amplimde  of  wavespeed  variability  is  diminished  below  200  km  depth  but 
significant  contrasts  remain.  At  200  km  depth  (Figure  6C)  the  images  reveal  again  a 
difference  between  the  wavespeeds  in  eastern  and  central  Australia.  The  eastern  edge 
of  the  high  wave  speeds  that  characterize  the  regions  of  Proterozoic  basement  seems  to 
coincide  with  the  southeastern  margin  of  the  Mt.  Isa  inlier  (20°S,  140°E)  and  the  western 
boundary  of  the  Eromanga  basin  and  extend  approximately  along  the  140°E  meridian 
down  the  eastern  side  of  the  Gawler  craton  (Figure  lA). 

At  depths  larger  than  300  km  (Figure  6D)  the  general  trends  in  the  anomalies  is 
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markedly  different  than  in  the  shallower  upper  mantle  (Figures  6A-C).  The  predominant 
orientation  of  the  gradients  is  no  longer  west-east,  although  even  at  these  large  depths 
the  strongest  low  wavespeed  anomalies  are  located  in  the  east,  that  is,  beneath  the 
Coral  Sea  and,  in  particular,  the  Tasman  Sea.  Between  300  and  400  km  depth  there 
are  prominent  high  wavespeeds  beneath  the  eastern  part  of  the  continent,  in  a  broad 
band  extending  from  20°S,  143°E  south  to  around  35°.  We  have  previously  noted  the 
high  velocities  beneath  the  New  England  fold  belt  near  30°S,  152°.  These  anomalies 
(discussed  in  more  detail  by  Zielhuis  &  Van  der  Hilst  1996)  are  robust  features  of  the 
solution  and  are  confirmed  by  a  simple  differential  technique  applied  to  a  small  number 
of  high-quality  Rayleigh  wave  data  (Passier  et  al.  1997).  The  high  wavespeeds  beneath 
the  western  parts  of  the  Musgrave  blocks  and  the  Officer  basin  may  be  real,  but  we  can 
not  be  confident  in  the  mapping  of  this  structure  untU  data  from  the  Skippy  arrays  in 
West  Australia  have  been  analyzed.  For  the  deeper  part  of  the  upper  mantle,  the  images 
display  several  ’streaks’  that  may  result  from  a  combination  of  uneven  sampling,  effects 
of  source  mislocation,  or  a  arise  from  the  neglect  of  inter-mode  coupling  for  the  higher 
mode  packet.  This  will  be  subject  to  further  study. 

Vertical  sections  through  the  three-dimensional  shear  wavespeed  model  for  the  upper 
mantle  (figure  7)  provide  further  insight  into  the  character  of  the  lithosphere  in  different 
regimes.  For  the  purpose  of  the  present  paper  we  loosely  associate  the  “lithosphere” 
with  a  zone  of  elevated  shear  velocity  compared  with  the  asthenosphere  beneath.  The 
locations  of  the  cross  sections  are  given  in  Figure  lA. 

A  section  from  the  Coral  Sea  across  the  eastern  margin  of  the  Australian  continent 
into  the  Proterozoic  shield  (figure  7A)  displays  the  large  contrast  in  lithosphere 
thickness  between  the  oceanic  and  continental  shield  regions.  Beneath  the  Coral  Sea 
a  negative  gradient  in  shear  velocity  starts  at  a  depth  of  about  50  km,  beneath  the 
coastal  region  of  easternmost  Australia  the  lithosphere  thickness  is  about  100  km,  and 
reduced  shear  wavespeeds  do  not  start  until  a  depth  of  about  2  50  km  beneath  the  central 
Proterozoic  cratons  (Figure  7A).  The  virtual  absence  of  a  high  wavespeed  lid  near  145°E 
coincides  with  the  region  of  Neogene  volcanism  in  the  Queensland  volcanic  province, 
which  may  suggest  that  part  of  the  continental  lithosphere  has  been  eroded  by  thermal 
processes. 

The  east-west  cross  section  at  24°S  (Figure  7B)  reveals  a  similar  contrast  in  lithosphere 
thickness  between  eastern  Australia  and  the  central  shields.  The  high  wavespeeds 
associated  with  the  seismological  lithosphere  are  detected  to  depths  exceeding  300 
km  beneath  West  Australia  (130°E).  The  image  clearly  reveals  the  low  velocity  zone 
extending  beneath  eastern  Australia  (east  of  140°E)  that  separates  thin  high  wavespeed 
lid  from  the  deeper  zone  of  high  wavespeed  anomalies  between  300  and  400  km  depth. 
This  cross  section  intersects  the  volcanic  province  at  approximately  150°E. 

The  seismic  structure  shallower  than  80  km  may  be  contaminated  by  variations  in 
crustal  structure  that  are  not  adequately  accounted  for  by  our  approach,  but  the  cross 
sections  clearly  reveal  a  pronounced  wavespeed  anomaly  in  the  region  associated  with 
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the  Alice  Springs  orogeny  and  the  eastern  part  of  the  Canning  basin  (125°-135°E).  At 
30°S  the  model  reveals  the  large  lateral  variations  in  wavespeed  near  the  Broken  Hill 
block,  i.e.  across  the  southern  part  of  the  Tasman  Line.  At  this  latitude  the  seismically 
defined  lithosphere  is  about  200  to  2  50  km  thick  beneath  most  of  the  Proterozoic  shields 
(between  125”  and  140°E),  and  possibly  even  thicker  beneath  the  Archaean  cratons 
further  west.  The  latter  observation  is,  however,  tentative  because  data  coverage  is  not 
sufficient  to  constrain  the  strucmres  in  the  southwestern  part  of  the  continent. 

The  map  views  (Figures  6A-C)  and  the  vertical  sections  (Figures  7A-C)  clearly  delineate 
the  transition  from  thin  hthosphere  beneath  the  Phanerozoic  of  eastern  Australia  to  the 
thick  Proterozoic  shields  of  central  Australia.  Locally,  in  particular  south  of  25”S.  the 
transition  is  complex  and  may  comprise  multiple  boundaries  and  basement  types. 

The  north-south  cross  sections  (Figures  7D-F)  further  illustrate  the  pronounced 
differences  in  character  of  the  high  wavespeed  lid  between  east  Australia  on  the  one 
hand  and  the  central  and  western  parts  of  the  continent  on  the  other.  The  thickest 
hd  is  located  beneath  the  western  part  of  the  Officer  basin  (at  about  25°S,  129”E) 
(Figure  7D),  that  is  near  the  boundary  of  the  Proterozoic  and  Archaean  shields.  The 
observation  of  such  a  thick  high  wavespeed  lid  in  the  western  part  of  the  continent, 
and  the  absence  of  a  strong  negative  velocity  gradient  in  shear  velocity  beneath  the 
lid  is  in  good  agreement  with  the  path  average  shear  wave  profile  by  Gaherty  &  Jordan 
(1995).  At  129°E  (Figure  7D)  the  sudden  increase  in  wavespeed  at  about  35”S  marks  the 
ocean-continent  transition  in  the  great  Australian  Bight,  and  the  deep  high  wavespeeds 
north  of  10”S  are  probably  related  to  subduction  beneath  eastern  Indonesia.  The  cross 
section  at  143°E  (Figure  7E)  displays  the  intermediate  thickness  of  the  lithosphere  in  the 
the  region  where  the  transition  from  the  coastal  region  to  the  central  shields  is  gradual. 
The  rapid  lateral  variations  near  30”S  coincide  with  the  boundary  between  the  Murray 
and  Eromanga  basins  (i.e.,  parts  of  the  Broken  Hill  block,  the  Lachlan  fold  belt,  and  the 
Darling  basin).  This  section  also  illustrates  the  deep  anomaly  near  22”S  in  Queensland. 
Figure  7F  displays  the  deep  low  wavespeed  anomaly  beneath  the  Tasman  sea  region,  the 
low  velocity  zone  beneath  easternmost  Australia,  and  the  deep  high  wavespeed  anomaly 
beneath  the  New  England  fold  belt. 

These  significant  results  have  been  derived  from  the  analysis  of  data  from  only  part 
of  the  Skippy  project.  More  detail  on  lithospheric  and  mantle  structure  beneath  the 
western  part  of  Australian  region  will  be  revealed  as  the  data  from  the  last  two  Skippy 
deployments  are  Incorporated  into  the  inversion  for  three-dimensional  structure. 


3.2  P  wave  tomography 

An  additional  source  of  information  on  three-dimensional  structure  comes  from  P  wave 
tomography  using  the  residuals  of  observed  arrival  times  of  P  phases  compared  with 
the  predictions  from  a  suitable  reference  model.  The  residuals  reflect  the  integrated 
influence  of  wavespeed  variations  in  the  Earth’s  interior  along  the  propagation  path 
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from  source  to  receiver.  With  sufficient  crossing  paths,  the  travel  time  information  can 
be  used  to  reconstruct  the  three-dimensional  variations  in  seismic  wavespeed. 

In  order  to  produce  a  homogeneous  set  of  travel  time  residuals  the  arrival  times 
determined  from  the  Skippy  records  are  incorporated  in  the  data  processing  scheme 
of  the  U.S.  Geological  Survey’s  National  Earthquake  Information  Center  (NEIC).  In  this 
way,  the  Skippy  data  help  constrain  the  hypocentres  of  earthquakes  in  the  Australasian 
region,  and  in  return  are  made  to  be  consistent  with  the  global  data  set  assembled  by 
Engdahl,  Van  der  Hilst  &  Buland  (1997). 

Travel  time  residuals  from  the  Skippy  project  were  first  used  by  Widiyantoro  &  Van 
der  Hilst  (1996)  in  their  tomographic  study  of  the  complex  subduction  beneath  the 
Indonesian  region.  Their  three-dimensional  mantle  model  was  generated  by  embedding 
a  high-resolution  representation  of  regional  structure  in  a  somewhat  lower  resolution 
cellular  model  for  the  rest  of  the  mantle.  By  this  means  contamination  of  the  regional 
structure  by  features  outside  the  region  of  interest  can  be  minimised.  The  inversions  for 
the  Indonesian  region  incorporated  phase  data  from  the  Skippy  records  in  (Queensland 
and  permanent  stations  in  Australia  and  reveal  the  thick  lithosphere  under  northern 
part  of  central  Australia. 

With  more  phase  picks  now  available  from  later  Skippy  deployments  (SKI,  SK2,  SK3, 
and  BAS),  in  addition  to  data  from  the  network  of  short-period  permanent  stations, 
the  P  wave  tomography  has  been  extended  to  cover  most  of  the  continent,  with 
resolution  in  the  eastern  part  significantly  better  than  in  the  west.  Figure  8  displays 
the  lateral  variation  in  P  wavespeed  relative  to  the  akl35  reference  model  at  a  depth 
of  approximately  150  km  beneath  the  Australian  region.  The  image  reveals  the  lateral 
contrast  between  high  and  low  wavespeeds  in  northern  Queensland  which  is  prominent 
in  the  shear  wave  structure  inferred  from  the  waveform  analysis  described  above.  The 
resolution  of  structure  in  oceanic  regions  is,  of  course,  poor  because  of  lack  of  station 
coverage. 

Much  of  the  information  used  to  construct  the  P  wave  velocity  model  comes  from 
arrivals  travelling  rather  steeply  in  the  mantle  so  that  vertical  resolution  is  limited. 
Nevertheless  the  strong  contrast  in  the  seismic  signature  of  the  lithosphere  between 
western  and  central  Australia  and  the  eastern  seaboard  are  clearly  revealed.  As  in 
the  shear  wave  images  the  zone  of  higher  wavespeeds  extends  further  east  than  the 
conventional  Tasman  Line  but  thicker  lithosphere  lies  mostly  to  the  west  of  140°E. 
There  are  some  significant  differences  in  the  P  and  S  wavespeed  models,  for  instance 
beneath  the  New  England  fold  belt.  These  differences  may  largely  be  due  to  variations 
in  data  coverage  and  resolution,  but  upon  further  analysis  they  may  also  yield  new 
information  about  the  physical  nature  of  the  anomalies,  as  for  example  whether  their 
origin  is  thermal  or  compositional. 
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3.3  Receiver  based  studies 

Along  with  their  use  in  tomographic  studies,  the  even  distribution  of  portable  broadband 
stations  across  Australia  can  also  be  exploited  to  provide  a  range  of  information  on  the 
structure  beneath  each  of  the  stations. 

3.3.1  Receiver  function  analysis 

For  each  of  the  broadband  stations  in  the  Skippy  deployments,  data  from  distant  seismic 
events  are  being  processed  to  extract  a  receiver  function  which  is  sensitive  to  the 
structure  in  the  crust  and  the  uppermost  mantle  and  particularly  the  character  of 
major  seismic  interfaces  such  as  the  Moho.  The  analysis  uses  the  P  wave  onset  and 
its  immediate  coda.  The  radial  and  vertical  components  of  the  record  share  the  same 
source  time  function,  so  that  the  dependence  on  the  radiation  from  the  source  can 
be  largely  eliminated  by  deconvolving  the  waveform  segment  of  the  radial  horizontal 
component,  lying  along  the  great-circle  path  from  the  source,  with  the  corresponding 
segment  of  the  vertical  component  (Langston  1979).  The  deconvolved  radial  receiver 
function  then  emphasises  the  influence  of  near  receiver  structure  and  can  be  inverted  in 
the  time  domain  for  a  1-D  shear  wave  velocity  model  of  the  crust  and  uppermost  mantle 
(see  e.g.  Ammon  et  al.  1990).  A  major  contribution  to  the  receiver  function  comes  from 
conversions  between  P  and  S  waves.  The  timing  and  amplimde  of  such  arrivals  provide 
constraints  on  the  properties  of  major  interfaces  such  as  the  crust-mantle  bovmdary 
and  internal  boimdaries  within  the  crust.  For  example,  a  sharp  crust/mantle  boundary 
produces  a  converted  phase  with  about  5  s  separation  from  P. 

The  inversion  of  the  receiver  functions  to  recover  crustal  and  uppermost  mantle 
structure  is  widely  recognised  to  be  sensitive  to  the  starting  model  if  a  conventional 
linearisation  scheme  is  employed  (see  e.g.  Ammon  et  al.  1990).  However,  such  difficulties 
can  be  overcome  by  employing  an  inversion  scheme  based  on  a  Genetic  Algorithm 
(Shibutani  et  al.  1996).  This  approach  makes  use  of  a  “cloud”  or  “population”  of  models 
to  minimise  the  dependence  on  a  starting  model;  a  set  of  biological  analogues  are  used 
to  produce  new  generations  of  models  from  previous  generations  with  preferential 
development  of  models  with  a  good  fit  between  observed  and  theoretical  receiver 
function's.  The  approach  provides  a  good  sampling  of  the  model  space  and  enables 
the  estimation  of  the  shear  wavespeed  distribution  in  the  crust  along  with  an  indication 
of  the  ratio  between  P  and  S  wavespeeds.  Many  models  with  an  acceptable  fit  to  data  are 
generated  during  the  inversion  and  a  stable  crustal  model  is  produced  by  employing 
a  weighted  average  of  the  best  1000  models  encountered  in  the  development  of  the 
genetic  algorithm.  The  weighting  is  based  on  the  inverse  of  the  misfit  for  each  model  so 
that  the  best  fitting  models  have  the  greatest  influence. 

At  each  station  a  stacked  receiver  fimction  was  constructed  from  teleseismic 
observations.  For  each  event,  the  receiver  function  was  low-pass  filtered  to  eliminate 
frequencies  above  1  Hz  in  order  to  minimize  the  influence  of  small  scale  heterogeneities. 
Subsequently,  the  receiver  functions  from  6-18  teleseismic  events  at  each  station 
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were  stacked  together  for  a  set  of  ranges  of  backazimuths.  The  weighting  used  in 
the  stacking  emphasizes  receiver  functions  with  higher  signal-to-noise  ratios  and  those 
whose  backazimuth  and  incident  angle  were  closer  to  the  mean  for  the  set  of  events. 

The  radial  component  of  the  stacked  receiver  function  was  then  inverted  for  a  1-D 
velocity  model  beneath  each  station.  In  the  Genetic  Algorithm  inversion,  the  crust  and 
uppermost  mantle  down  to  60  km  was  modelled  with  six  major  layers:  a  sediment 
layer,  basement  layer,  upper  crust,  middle  crust,  lower  crust  and  uppermost  mantle. 
The  model  parameters  in  each  layer  are  the  thickness,  the  S  wave  velocity  at  the  upper 
boundary,  the  S  wave  velocity  at  the  lower  boundary  and  the  velocity  ratio  between  P  and 
S  waves  (VpIVs).  The  S  wave  velocity  for  each  layer  is  constructed  by  linearly  connecting 
the  values  at  the  upper  and  lower  boundaries,  to  give  a  sequence  of  constant  velocity 
gradient  segments  separated  by  velocity  discontinuities. 

The  results  of  the  receiver  function  inversion  for  the  stations  SB06  and  SAG  7  are 
shown  in  figure  9.  All  10,000  models  searched  in  the  Genetic  Algorithm  inversion  for 
5  velocity  are  shown  as  the  light  gray  shaded  area  and  superimposed  on  this  the  best 
1,000  models  are  shown  with  darker  gray  tones,  where  the  darkness  is  logarithmically 
proportional  to  the  number  of  the  models.  The  best  fitting  model  for  each  site  is  shown 
as  a  black  line.  However  a  more  useful  and  stable  result  is  provided  by  the  averaged 
model  generated  by  weighting  the  best  1,000  models  by  the  inverse  of  their  misfit 
values.  This  averaged  model  is  shown  in  white  for  S  velocity  and  in  a  solid  line  for  the 
VpjVs  ratio.  The  lower  panel  in  figure  9  compares  the  waveforms  of  the  stacked  receiver 
function  to  synthetics  calculated  for  the  averaged  models.  The  fit  to  the  major  phases  is 
good,  and  this  will  be  true  of  most  of  the  best  1, 000  models.  The  zone  of  darker  gray  in 
the  upper  part  of  figure  9  can  therefore  be  thought  of  as  an  indicator  of  the  constraint 
on  the  crust  and  upper  mantle  structure.  Both  SB06  and  SA07  have  surficial  sediment 
with  very  low  velocities  so  that  P-to-S  converted  phases  and  reverberations  (1  -  3  s) 
originating  in  the  sediment  layer  are  larger  than  the  direct  P  phase  which  arrives  at  zero 
time. 

The  analysis  for  crustal  structure  has  been  completed  for  all  the  Skippy  stations  in  the 
SKI  and  SK2  arrays  in  eastern  Austraha  and  is  in  progress  for  the  permanent  stations 
and  the  later  Skippy  arrays.  The  parametrisation  of  the  model  at  each  station  is  via  a 
sequence  of  velocity  gradients  and  discontinuities;  sensitivity  analysis  indicates  that  the 
dependence  on  phase  conversions  enhances  the  resolution  of  boundaries  at  the  expense 
of  the  gradients.  Nevertheless,  all  the  models  have  been  derived  by  the  same  procedure 
and  we  can  make  direct  comparisons  and  extract  a  wide  range  of  information  on  crustal 
structure  and  uppermost  mantle  structure  across  eastern  Australia.  In  the  few  places 
where  direct  comparisons  can  be  made,  e.g.  ZB12,  there  is  a  very  good  concordance 
between  the  character  of  the  S  velocity  structure  from  the  receiver  function  inversion 
and  the  previous  P  velocity  model  from  refraction  studies. 

Figures  10,  11  and  12  summarise  the  properties  of  the  crustal  models  at  each  of 
the  Skippy  stations  in  eastern  Australia  in  terms  of  the  major  subdivisions  of  the  crust: 
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the  upper  crust  (fig  10),  middle/lower  crust  (fig  11)  and  the  crust-mantle  boundary 
and  uppermost  mantle  velocities  (fig  12).  In  each  map  we  indicate  the  thickness  of 
the  layer,  the  corresponding  crustal  velocities  and  the  character  of  the  boundary.  A 
sharp  boundary  is  associated  with  a  clear  converted  phase  in  the  receiver  function 
waveform  whereas  the  expression  of  a  transitional  zone  is  more  subtle.  The  nature 
of  the  boundaries  are  classified  into  the  categories:  sharp  (<  1  km),  thin  (<  4  km), 
TRANSITIONAL  Or  INTERMEDIATE  (<  10  km),  and  BROAD  (>  10  km).  The  depth  of  the  crustal 
boundaries  was  estimated  at  each  station;  for  a  transitional  zone  the  lower  boundary 
was  selected. 

The  Information  from  the  receiver  function  inversion  provides  a  major  supplement 
to  the  previous  results  from  isolated  refraction  and  reflection  experiments  and  in  the 
crustal  thickness  map  in  figure  12  the  earlier  information  (from  Collins  1991,  indicated 
by  crosses)  has  been  combined  with  the  present  results.  The  crust-mantle  boundary  is 
deep  (38  -  44  km)  and  mostly  transitional  in  character  along  the  axis  of  the  fold  belt 
zone  in  the  east  (from  stations:  ZB12,  SB09  to  SA05,  SA03).  A  relatively  sharp  Moho  is 
found  at  a  shallower  depth  (30  -  36  km)  at  the  western  edge  of  the  study  area  close  to 
the  boundary  between  Phanerozoic  and  Precambrian  exposure. 

A  major  advantage  of  the  receiver  function  approach  is  that  it  provides  good 
constraints  on  shear  velocities  in  the  crust  and  uppermost  mantle,  which  are  not  well 
resolved  by  the  partitioned  waveform  inversion  of  the  S  wave  and  surface  wave  portions 
of  the  seismogram.  The  3-D  models  are  most  reliable  from  60  km  down,  whereas  the 
receiver  function  analysis  is  most  effective  above  60  km.  There  is  a  good  agreement 
between  the  results  from  the  waveform  tomography  and  the  receiver  function  analysis 
at  60  km  depth. 


3.3.2  SKS  splitting 

In  addition  to  mapping  isotropic  variations  in  shear  wavespeed,  the  Skippy  project  has 
also  enabled  us  to  investigate  continental  azimuthal  anisotropy  on  an  unparalleled  scale. 
From  observations  of  shear  wave  splitting  of  the  seismic  core  phase  SKS,  the  angle  of  fast 
polarisation  direction  (<^)  and  the  time  difference  between  fast  and  slow  SKS  waves  (dt) 
has  been  extracted  for  many  of  the  Skippy  stations  and  a  number  of  permanent  stations, 
see  Clitheroe  &  Van  der  Hilst  (1997,  this  issue).  The  broad-band  data  demonstrate  that 
SKS  splitting  is  manifested  at  different  frequencies.  At  a  frequency  of  about  1  Hz, 
dt  is  0.3— 0.6  s,  and  <f)  parallels  pre-existing  crustal  fabric.  At  lower  frequencies  </> 
varies  across  the  continent:  0  is  inconsistent  with  present-day  plate  motion  direction, 
but  can  probably  be  explained  by  predominant  mineral  orientation  in  the  sub-crustal 
lithosphere.  The  resemblance  between  the  location  of  the  assumed  eastern  boundary  of 
the  central  cratons  and  the  spatial  pattern  of  (/)  in  the  east  of  Australia  may  suggest  that 
shear  wave  splitting  carries  information  about  either  the  deep  structure  of  the  mountain 
belts  formed  during  the  Palaeozoic  accretion  of  terranes  onto  the  central  cratons  or 
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about  astheno spheric  flow  around  the  undulating  eastern  edge  of  the  deep  cratonic 
keel. 


4  CONCLUSION 

Over  the  last  decade  there  has  been  a  significant  increase  in  knowledge  of  the  P  and 

5  wave  velocities  in  the  mantle,  particularly  beneath  northern  Australia.  The  current 
generation  of  three-dimensional  studies  based  on  the  use  of  portable  seismometers  have 
the  potential  to  dramatically  increase  the  level  of  understanding  of  mantle  structure 
beneath  the  Australian  region  and  in  particular  the  contrast  between  eastern  and 
western  Australia. 
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Figure  1.  (A)  Map  of  the  crustal  elements  and  tectonic  units  referred  to  in  the  text,  together 
with  with  the  lines  of  the  sections  in  figure  7.  [Key:  Ar  -  Arunta  Block,  Am  -  Amadeus  Basin,  BH  - 
Broken  HiU  Block,  Ca  -  Canning  Basin,  Da  -  Darling  Basin,  Er  -  Eromanga  Basin,  Eu  -  Eucla  Basin, 
Ga  -  Gawler  Block,  HC  -  HaUs  Creek  Belt,  Ki  -  Kimberley  Block,  La  -  Lachlan  Fold  Belt,  MI  -  Mt  Isa 
Block,  Mu  -  Musgrave  Block,  NE  -  New  England  Fold  Belt,  Or  -  Ord  Basin,  Of  -  Officer  Basin].  (B) 
Configuration  of  portable  seismic  recording  stations  1985-1996.  Short-period  stations  are  indi¬ 
cated  by  open  triangles.  Portable  broadband  stations  are  marked  by  sohd  symbols.  Permanent 
stations  with  high  fideUty  recording  are  indicated  by  a  double  circle  and  station  name.  The  ap¬ 
proximate  location  of  boundary  between  Precambrian  outcrop  in  western  and  central  Australia 
and  the  Phanerozoic  east  is  indicated  by  a  shaded  line. 

Figure  2.  A  composite  record  section  of  unfiltered  broadband  seismograms  for  five  New  Guinea 
events  recorded  at  Skippy  stations  in  Queensland.  The  section  is  constructed  from  the  vertical 
component  for  each  path  and  covers  the  P  and  SV  waves  returned  from  the  upper  mantle.  The 
travel  time  curves  are  for  an  event  at  25  km  depth  in  the  akl35  reference  model  (Kennett, 
Engdahl  &  Buland  1995) 


Figure  3.  Enlargement  of  the  S  wave  arrivals  for  composite  record  sections  of  unfiltered  broad¬ 
band  seismograms  for  five  New  Guinea  events  recorded  at  Skippy  stations  in  Queensland.  Sec¬ 
tions  are  shown  for  both  the  radial  (R)  and  tangential  components  (T)  constructed  for  each  path. 
The  travel  time  curves  are  for  an  event  at  25  km  depth  in  the  akl35  reference  model  (Kermett, 
Engdahl  &  Buland  1995) 

Figure  4.  P  and  S  velocity  and  attenuation  structure  for  the  upper  mantle  beneath  northern 
Australia  determined  from  a  combination  of  short-period  and  broadband  observations 

Figure  5.  Wavepath  coverage  available  for  the  Skippy  experiment  in  eastern  and  central  Australia, 
(a)  wave  paths  exploited  by  Zielhuis  &  Van  der  ffilst  (1996);  (b)  wave  path  covergae  use  in  this 
study. 


Figure  6.  Three-dimensional  shear  wavespeed  model  derived  from  partitioned  waveform  inver¬ 
sion  from  the  Skippy  experiment  in  eastern  and  central  Australia.  Map  views  at  (a)  80  km,  (b)  140 
km,  (c)  200  km,  and  (d)  300  km  depth.  Based  on  the  path  coverage  (Figure  5B)  the  wavespeed 
anomalies  are  reUably  images  east  of  about  125°E. 

Figure  7.  Vertical  cross-sections  through  the  three-dimensional  shear  wavespeed  derived  from 
partitioned  waveform  inversion  from  the  Skippy  experiment  in  eastern  and  central  Australia. 
Cross  sections  at  (a)  lO'S,  (b),  24“S,  and  (c)  30“S,  and  north-south  sections  at  (d)  129'-E,  (e) 
143'’E,  and  (f)  ISl’E 

Figure  8.  Cross-section  through  the  3D  model  of  the  P  wave  velocities  in  the  Australian  region 
for  the  depth  interval  from  160-220  km 
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Figure  9.  Genetic  algorithm  inversion  of  a  receiver  function  to  determine  S  wave  structure  and 
the  Vp/Vs  ratio.  All  10,000  models  searched  in  the  GA  inversion  are  shown  as  the  light  gray 
shaded  area.  The  best  1,000  models  are  shown  with  darker  gray  shaded  area.  The  darkness 
is  logarithmically  proportional  to  the  number  of  the  models  as  shown  by  the  scale  bar.  The 
best  model  and  the  averaged  model  are  shown  by  the  black  solid  line  and  the  white  solid  line, 
respectively.  For  the  Vp/Vs  ratio,  the  the  solid  line  indicates  the  averaged  model. 

Figure  10.  Summary  of  velocity  and  thickness  information  for  the  upper  crust  in  eastern  Aus¬ 
tralia 

Figure  11.  Summary  of  velocity  and  thickness  information  for  the  middle  and  lower  crust  in 
eastern  Australia 

Figure  12.  The  depth  of  the  crust-mantle  boundary  beneath  eastern  Australia,  and  the  S  wave 
velocity  and  Vp/Vs  ratio  at  the  top  of  the  mantle.  The  nature  of  the  crust-mantle  boundary  is 
classified  into  four  categories  :  sharp  (<  1  km),  thin  (<  4  km),  intermediate  (<  10  km),  and 
BROAD  (>  10  km).  The  crosses  indicate  the  locations  for  which  the  crustal  structure  has  been 
obtained  previously  from  refraction  surveys  [Collins,  1991].  These  data  were  incorporated  into 
the  contouring. 
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Part  11:  Upper  mantle  heterogeneities  in  the  Indian  Ocean  from  waveform  inversion. 
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Abstract.  A  waveform  inversion  method  is  applied  to 
156  Love  and  Rayleigh  wave  seismograms  to  build  up 
a  3-D  model  of  the  shear  velocity  in  the  upper  mantle 
beneath  the  Indian  Ocean.  The  first  step  of  the  method 
consists  in  finding,  for  each  path,  a  radially  strati¬ 
fied  upper  mantle  model  compatible  with  the  Love  and 
Rayleigh  wave  seismograms  relative  to  that  path.  The 
fundamental  mode  and  few  higher  modes  are  modelled 
in  the  waveform  inversion.  In  a  second  step,  the  models 
related  to  the  different  paths  are  used  in  a  tomographic 
inversion  to  map  the  3-D  upper  mantle  structure.  The 
3-D  velocity  model  has  a  lateral  resolution  of  1000  km. 
No  significant  velocity  anomalies  are  found  below  300 
km,  although  the  resolution  is  still  good.  Continental 
roots  are  confined  to  the  upper  300  km  and  low  ve¬ 
locities  below  mid-oceanic  ridges  to  the  upper  250  km, 
depending  on  their  spreading  rates.  At  shallow  depths 
(<  80  km),  we  find  a  positive  velocity  anomaly  beneath 
the  West  Indian  ridge  near  the  Rodriguez  triple  junc¬ 
tion,  where  gravimetric  and  bathymetric  data  indicate 
a  reduced  volcanic  activity.  At  greater  depth  (around 
200  km),  a  low  velocity  signature  is  found  beneath  the 
hot-spot  of  Reunion-Mauritius  islands  and  the  Central 
Indian  ridge.  It  could  reflect  a  real  connection  between 
the  two  structures. 


1.  Introduction 

During  the  past  15  years,  long  period  surface  waves 
have  been  commonly  used  to  image  the  structure  of  the 
upper  mantle.  Due  to  the  increasing  amount  of  digi¬ 
tal  data,  the  lateral  dimensions  of  structural  features 
that  can  be  resolved  in  global  S  velocity  models  have 
progressively  decreased  from  5000  km  [  Woodhouse  and 
Dziewonski,  1984]  to  around  1000  km  [Zhang  and  Tan- 
imoto,  1993].  In  regional  studies,  lateral  resolution  can 
now  reach  locally  250  km  in  regions  where  a  dense  cover¬ 
age  is  available  [Zielhuis  and  van  der  Hilst,  1996].  In  the 
Indian  Ocean,  the  previous  regional  studies  [Montagner^ 
1986a;b;  Roult  et  al,,  1987;  Montagner  and  Johert,  1988] 
have  a  lateral  resolution  limited  to  2000  km,  due  to 
the  poor  distribution  of  stations.  These  latter  studies, 
based  on  the  analysis  of  fundamental  mode  dispersion 
curves,  show  a  good  correlation  of  the  deep  structure 
with  surface  tectonics  down  to  about  100  km  depth,  ex¬ 
cept  in  the  Central  part  of  the  Southeast  Indian  Ridge. 
In  addition,  Montagner  [1986a;b]  and  Montagner  and 
Johert  [1988]  find  a  slow  structure  shifted  eastward  of 
the  Central  Indian  ridge  below  100  km  depth. 

In  this  paper,  we  present  a  3-D  S-velocity  model  for 
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the  upper  mantle  in  the  Indian  Ocean  with  improved 
lateral  and  depth  resolution.  The  model  is  obtained 
from  a  larger  dataset  than  in  previous  studies  and  a 
different  methodological  approach.  New  data  recorded 
on  Geoscope  stations  DRV,  CRZ  and  HYB  allow  us  to 
improve  the  lateral  resolution  to  1000  km  with  a  more 
homogeneous  path  coverage.  Instead  of  analysing  dis¬ 
persion  curves,  we  use  a  waveform  inversion  method 
designed  to  match  directly  surface-wave  seismograms, 
possibly  including  several  higher  modes.  It  is  based  on 
the  waveform  inversion  technique  of  Cara  and  Leveque 
[1987]  adapted  to  process  simultaneously  a  set  of  seis¬ 
mograms  with  close  epicenters,  recorded  at  a  single  sta¬ 
tion.  In  this  approach  called  ’’multi-seismogram  tech¬ 
nique”,  the  epicenters  of  a  given  cluster  must  be  close 
enough  so  that  we  can  consider  the  waves  travel  the 
same  path  and  can  be  explained  with  a  single  model. 
The  inverted  model,  compatible  with  all  the  observed 
seismograms  available  for  that  path,  can  be  interpreted 
as  the  average  Earth’s  structure  beneath  the  path.  A 
regionalization  of  the  models  obtained  for  the  different 
paths  leads  to  tomographic  images  of  the  upper  mantle. 


2.  Data  Selection 


We  use  Love  and  Rayleigh  wave  long-period  seismo¬ 
grams  at  11  stations  in  the  Indian  Ocean,  for  earth¬ 
quakes  occurring  from  1988  to  1994.  The  selected  events 
have  magnitudes  ranging  from  4.8  to  7.2  Msz  and  the 
corresponding  seismograms  are  inverted  using  observ¬ 
ables  covering  the  20-300  s  period  range.  Most  of  these 
earthquakes  are  shallow  events  located  on  mid-oceanic 
ridges  of  the  Indian  Ocean,  The  few  deep  earthquakes 
we  use  are  located  in  the  Indonesian  subduction  zone. 
We  have  discarded  seismograms  with  a  poor  signal-to- 
noise  ratio  or  for  which  the  initial  phase  at  the  source 
was  not  stable  with  respect  to  small  perturbations  of 
the  path  azimuth.  Finally,  we  keep  only  the  paths  for 
which  at  least  one  Rayleigh  wave  and  one  Love  wave 
seismograms  can  be  inverted  simultaneously.  For  some 
paths,  up  to  4  seismograms  (2R4-2L)  have  been  inverted 
simultaneously.  The  whole  set  of  data  consists  of  156 
seismograms  related  to  71  different  paths.  For  11  of 
these  paths,  higher  modes  are  taken  into  account  up  to 
the  fifth  overtone  for  Rayleigh  waves.  Love  wave  higher 
modes  are  included  in  the  inversion  for  53  paths.  The 
path  coverage  is  shown  on  Fig.  la. 


Fig.  la 


3,  Method 

3.1.  The  Waveform  Fitting 

One  important  problem  in  a  waveform  inversion  tech¬ 
nique  is  the  highly  non-linear  dependence  of  the  recorded 
signal  on  the  model  parameters.  The  originality  of  the 
method  developped  by  Cara  and  Leveque  [1987],  com¬ 
pared  with  other  waveform  inversion  techniques  [e.g. 
Nolet,  1990],  is  the  introduction  of  secondary  observ¬ 
ables,  built  up  from  the  seismograms,  having  only  a 


slightly  non-linear  dependence  upon  the  model  parame¬ 
ters,  allowing  for  an  inversion  with  a  classical  non-linear 
scheme.  In  the  multi- seismogram  approach  used  here, 
the  secondary  observables  are  computed  for  all  seismo¬ 
grams  related  to  a  same  path,  and  then  inverted  using 
the  same  scheme  as  in  Cara  and  Leveque  [1987]. 

Using  several  seismograms  simultaneously  in  one  in¬ 
version  reduces  the  influence  of  errors  on  fixed  parame¬ 
ters,  such  as  focal  mechanism  or  focal  depth,  and  allows 
to  use  simultaneously  data  containing  non-redundant 
information  on  the  mantle  structure.  For  example,  we 
can  invert  simultaneously  data  from  shallow  and  deep 
earthquakes  as  in  Stutzmann  and  Montagner  [1993],  or 
data  from  Love  and  Rayleigh  waves. 

In  most  cases,  Love  and  Rayleigh  waves  are  not  com¬ 
patible  with  a  single  isotropic  structure.  We  adress 
this  Love/Rayleigh  discrepancy  by  inverting  for  the  S- 
wave  anisotropic  parameter  ^  as  defined  in  Takeuchi 
and  Saito  [1972]  for  a  transversely  isotropic  medium. 
We  also  invert  for  the  shear  wave  velocity  the  at¬ 
tenuation  factor  and  for  Mo,  the  scalar  seismic 
moment  of  the  source,  to  ensure  an  adequate  scaling  of 
the  signal.  In  this  paper,  we  present  and  discuss  results 
on  pv  only.  Significant  anisotropy  is  found  in  the  up¬ 
per  250-300  km,  but  the  trade-off  between  Pv  and  ^  is 
found  very  small.  Results  on  anisotropy  and  will 
be  discussed  in  a  separate  paper. 

The  a  priori  information  used  in  this  study  is  similar 
to  the  one  used  in  Leveque  et  at.  [1991]:  the  a  pri¬ 
ori  Py  model  in  particular  is  a  smoothly  depth- varying 
isotropic  model  in  the  mantle  with  a  crustal  structure 
adapted  for  each  path.  The  a  priori  standard-deviations 
are  fixed  to  0.05  km/s  for  pv  and  0.03  for  ^  according  to 
the  expected  variation  range  for  these  parameters  in  the 
Indian  Ocean.  The  vertical  correlation  length  is  fixed 
to  50  km  and  the  uncertainty  on  the  data  to  10  %  for 
amplitude  data,  and  5  %  for  phase  data. 

3.2.  Regionalization  of  Shear  Velocities 

The  average  Pv  models  are  obtained  for  each  path 
shown  in  Fig. la.  A  regionalization  of  these  models 
is  performed  using  the  continuous  regionalization  algo¬ 
rithm  of  Montagner  [1986a].  To  constrain  the  lateral 
smoothness  of  the  model  and  according  to  our  path  cov¬ 
erage,  we  use  gaussian  functions  with  length  of  1000  km 
as  correlation  functions.  The  a  priori  standard  devia¬ 
tion  which  controls  the  amplitude  of  heterogeneities  in 
the  inverted  model  is  set  to  0.10  km/s.  This  value  al¬ 
lows  velocity  contrasts  in  the  final  model  as  large  as 
15  %  ,  in  agreement  with  velocity  contrasts  observed 
in  previous  oceanic  studies  [e.g.  Montagner,  1986a;b  or 
Nishimura  and  Forsyth  1989]. 

The  S- velocity  map  at  50  km  depth  and  an  East- West 
cross  section  at  latitude  20° S  are  shown  on  Fig.l,  to¬ 
gether  with  the  corresponding  a  posteriori  error  maps 
(Fig.la.b).  Dark  blue  areas  indicate  an  a  posteriori  er¬ 
ror  close  to  the  a  priori  error  (0.1  km/s),  indicating  a 
total  lack  of  resolution.  This  unresolved  area  is  limited 


to  the  outer  zone  of  the  map.  In  most  part  of  the  Indian 
Ocean,  and  at  least  down  to  300  km  depth,  the  error  is 
less  than  0.05  km/s,  indicating  a  very  good  resolution. 
Below  300  km  depth  (Fig. lb),  we  still  observe  locally  an 
error  smaller  than  0.05  km/s  in  two  regions,  one  located 
below  the  Central  Indian  ridge  (60°jF),  and  the  other 
one  around  llO^F?.  They  correspond  to  areas  where  the 
paths  associated  to  deep  events  in  the  eastern  part  of 
the  indonesian  subduction  zone  cross  the  latitude  20^5. 
This  is  a  direct  consequence  of  the  increased  resolution 
for  these  paths,  as  compared  with  paths  where  no  higher 
mode  data  were  available. 

At  50  km  depth  (Fig.lc),  S-velo cities  are  strongly 
correlated  to  surface  tectonics.  Low  velocities  are  as¬ 
sociated  with  mid- oceanic  ridges.  Old  oceanic  basins 
display  high  velocities.  Two  unexpected  features  how¬ 
ever  show  up  in  the  well  resolved  area:  a  positive 
anomaly  beneath  the  West  Indian  ridge,  near  the  Ro¬ 
driguez  triple  junction,  and  high  velocities  beneath  the 
Carlsberg  ridge.  On  the  cross-section  at  latitude  20^19 
(Fig. Id),  fast  S-velocities  are  found  down  to  250  km  be¬ 
low  Africa  and  down  to  300  km  below  Australia.  The 
low  velocity  signature  of  the  Central  Indian  ridge  does 
not  extend  beyond  250  km  depth.  This  is  probably  a 
reliable  feature  in  the  model  since  the  resolution  is  good 
in  this  area  (Fig. la).  Below  300  km,  the  amplitude  of 
the  anomalies  remains  small  everywhere.  Note  also  that 
the  low  velocity  beneath  the  island  of  Madagascar  can 
be  an  artefact  due  to  a  poor  knowledge  of  the  crustal 
structure  in  this  region;  this  negative  anomaly  does  not 
extend  deeper  than  50  km. 

Figure  2  displays  the  velocity  maps  at  increasing 
depths,  showing  that  features  observed  on  the  20°i9 
cross-section  are  more  general.  The  correlation  with 
surface  tectonics  is  still  very  clear  at  88  and  140  km 
depth.  At  88  km,  a  positive  velocity  contrast  is  present 
beneath  the  West  Indian  ridge  near  the  Rodriguez  triple 
junction.  At  140  km  depth,  this  positive  contrast  dis¬ 
appears  and  the  slow  signature  of  the  West  Indian  ridge 
itself  vanishes  while  it  remains  beneath  the  two  other 
ridges.  At  190  km  depth,  low  velocities  beneath  mid- 
oceanic  ridges  and  the  Indonesian  arc  nearly  disappear. 
At  this  depth,  the  largest  negative  anomaly  is  centered 
on  the  Central  Indian  ridge,  close  to  Mauritius  and  La 
Reunion  islands.  High  velocities  are  still  present  be¬ 
neath  continents,  and  a  north-south  positive  anomaly 
connects  the  Australian  continent  to  Antarctica.  At  268 
km,  the  low  velocity  signature  of  the  ridge  completely 
vanishes,  but  the  high  velocity  pattern  is  similar  to  that 
found  at  190  km  depth,  although  if  attenuated. 


Figure  2 


4,  Discussion 

The  depth  extent  of  the  seismic  anomalies  we  have 
found  (Fig. Id, 2c)  is  slightly  shallower  than  in  most 
other  surface  wave  studies  of  the  Indian  Ocean  [e.g. 
Montagner  and  Tanimoto,  1991].  These  observations 
agree  with  a  more  recent  global  S-velocity  model  [J. 
Woodhouse  and  J.  Trampert,  submitted  to  Earth  Planet. 


Sci.  Lett.,  1996]  and  support  the  idea  developed  in 
Ricard  et  al.  [1996]  that  most  seismological  observa¬ 
tions  are  compatible  with  the  3-SMAC  ’’geodynami- 
cal”  model  where  upper  mantle  anomalies  are  no  deeper 
than  300  km.  However,  in  our  model,  the  low- velocity 
signature  beneath  mid- oceanic  ridges  is  apparent  down 
to  200  km  in  the  most  part  of  the  Indian  Ocean  (Fig. 
2c, d),  a  greater  value  than  in  3-SMAC  where  ridge 
anomalies  are  present  in  the  first  100  km  only. 

Despite  the  enhanced  lateral  resolution  of  our  study, 
no  eastward  shift  is  observed  at  any  depth  beneath 
the  Central  Indian  ridge:  we  find  low  velocities  in  the 
Central  Indian  basin  at  190  km  depth  as  in  Montagu 
ner  [1986a;b]  or  Montagner  and  Jobert  [1988]  but  the 
strongest  velocity  anomaly  is  well  centered  beneath  the 
Central  Indian  ridge  (Fig.2).  Another  feature,  more 
classical  [e.g.  Montagner  and  TanimotOj  1991],  of  our 
model  is  the  correlation  between  the  low  velocity  signa¬ 
ture  of  oceanic  ridges  and  their  spreading  rates. 

At  shallow  depths,  we  find  a  positive  anomaly  that 
does  not  appear  in  previous  3D  S- velocity  models  be¬ 
neath  the  West  Indian  ridge  near  the  Rodriguez  triple 
junction,  suggesting  a  relatively  cold  mantle  in  that  part 
of  the  ridge  rFig.lc,2a).  This  can  be  related  to  recent 
bathymetric  [Mendel  and  Sauter,  1996]  and  gravimet¬ 
ric  [Rommevaux  et  alj  submitted  to  Marine  Geophys. 
Res.,  1996]  results  that  conclude  to  a  weak  volcanic 
activity  between  the  fracture  zone  of  Melville  (around 
60^E)  and  the  Rodriguez  triple  junction.  It  could  also 
help  explaining  the  high  phase  velocities  obtained  at 
this  location  by  Roult  et  al.  [1987]  for  Rayleigh  waves 
at  periods  61s  and  102s.  Another  interesting  new  fea¬ 
ture  is  the  deep  low  velocity  anomaly  that  persists  to 
190  km  depth  beneath  the  Central  Indian  ridge,  and 
spreads  toward  the  hot-spot  of  Reunion-Mauritius  is¬ 
lands.  This  westward  spreading  is  also  apparent  on  the 
cross-section  (Fig. Id)  and  at  shallower  depth  (88  and 
140  km,  Fig.2a,b)  beneath  the  Mascareignes  plateau. 
Due  to  a  lack  of  lateral  resolution,  this  low  velocity 
anomaly  could  be  the  coalescence  of  two  distinct  anoma¬ 
lies  beneath  the  hot-spot  of  Reunion-Mauritius  islands 
and  beneath  the  Central  Indian  ridge.  Another  possi¬ 
ble  interpretation  is  to  associate  this  low  velocity  signa¬ 
ture  with  the  hot-spot  trace  between  the  Central  Indian 
ridge  and  Mauritius  island.  The  well-known  AustreJian- 
Ant arctic  discordance  appears  clearly  on  our  maps  at  all 
depths  as  a  structure  abnormally  fast  for  a  ridge  (Fig.2). 
This  result  is  in  agreement  with  those  of  Forsyth  et  al 
[1987],  who  found  high  S  velocities  down  to  150  km 
depth,  and  Kuo  et  al.  [1996]  who  favour  a  deep  as- 
thenospheric  origin  of  the  AAD  from  considerations  on 
the  geoid.  The  oceanic  plateaus,  which  are  one  of  the 
striking  characteristics  of  the  Indian  Ocean,  do  not  sys¬ 
tematically  correspond  to  clear  velocity  anomalies.  This 
may  be  due  to  the  fact  that  some  of  these  structures 
(e.g.  Ninety-east  and  Broken  ridges)  are  much  narrower 
than  the  lateral  resolution  or,  where  the  plateau  struc¬ 
ture  is  larger  (Kerguelen),  to  the  poor  path  coverage 
near  the  border  of  our  maps. 
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Figure  1.  a)  A  posteriori  error  map  (in  km/s)  at  50 
km  depth  and  path  coverage,  b)  vertical  East- West 
cross-section  at  20^5  in  the  3-D  error  model,  c)  Map 
of  lateral  perturbations  in  percent  (with  respect  to  the 
value  Vref  indicated  on  the  map)  for  the  S-velocity  at  50 
km  depth,  d)  vertical  East-West  cross-section  at  20"^^ 
in  the  3-D  S-velocity  model.  The  oceanic  ridges  of  the 
Indian  Ocean  are  plotted  in  green  on  the  maps. 
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Figure  2.  Maps  of  lateral  perturbations  for  the  S- 
velocity  in  the  Indian  Ocean  at  increasing  depths.  The 
variations  are  given  in  percents  with  respect  to  the  value 
Vref  indicated  on  each  map  a)  depth  88  km.  b)  depth 
140  km.  c)  depth  190  km.  d)  depth  268  km. 
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Part  III:  Guided  waves  in  3-dimensional  structures 
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differential  equations  can  be  constructed,  for  the  components  of  displacement  and  the 

traction  in  the  preferred  direction.  Derivatives  with  respect  to  the  preferred  coordinate  Tij  =  CijkidkUi,  (2.1.4) 

are  restricted  to  first  order  only  and  can  be  separated  from  the  remaining  terms  in 

can  be  written  in  terms  of  traction  vectors  T,  with  components  (T,),  =  T,;,  in  the  form 

which  derivatives  with  respect  to  the  other  spatial  coordinates  and  time  appear.  This 

approach  is  well  suited  to  models  in  which  the  structure  is  two-dimensional.  Kennett  t,  =  C,/3;U.  (2.1. S) 


Guided  waves  in  3-dimensional  structures  11  12  B.LM  Kennett 


g  £  g 

^  Ti  B 


o  ■£ 

(sj  ^ 


o  *£ 


P<  •n  ^ 

^  S  § 


C  *-  a 

•c  5 

c  fa 
aj  5  fa 


;s  (V 

a  ^ 

B  -r 

.5  oj  H 

«  .a  ^ 

c  «  a 

eo  Rj  > 


I  £i  & 

H  «  aj 


We  consider  an  expansion  (2.2.1)  in  terms  of  the  modes  of  a  fixed  reference  structure  In  this  case,  the  modal  field  used  in  the  representation  (2.2.1)  is  chosen  to  match 

with  propagation  properties  described  by  the  differential  operators  The  same  set  the  properties  of  the  structure  along  a  vertical  profile  beneath  each  surface  point, 

of  eigenmodes  u®  is  used  throughout  and  so  there  will  be  no  contributions  from  the  The  eigenfunctions  u5(Xj.,X3)  will  therefore  vary  with  horizontal  position,  but  their 

horizontal  derivatives  of  the  modal  eigenfunctions  (i.c.  J  =  0).  properties  will  be  such  as  to  annihilate  any  contribution  from  the  integrals  over  the 

The  horizontal  evolution  equations  for  the  modal  coefficients  are  then  terms  involving  the  differential  operators. 


representing  an  orthogonality  relation  between  different  modes  for  an  anisotropic  and 
attenuative  halfspace.  For  a  perfectly  elastic  medium  we  can  recover  the  form  used 
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Abstract 

Seismic  velocity  heterogeneity  in  the  Earth’s  mantle  is  strongly  concentrated  near  its  top.  The 
shallow  heterogeneity  of  the  mantle  correlates  strongly  with  surface  tectonics.  We  use  these 
observations  as  constraints  of  a  tomographic  experiment  aimed  at  building  a  regionalized 
upper-mantle  (RUM)  reference  model.  We  use  a  select  set  of  teleseismic  travel  times  to  minimize 
the  mapping  of  mislocation  into  structure.  The  data  selection  emphasizes  the  robustness  of 
individual  picks.  The  form  of  the  RUM  model  is  a  set  of  velocity  profiles  as  functions  of  depth 
through  the  upper  mantle  for  each  of  the  different  tectonic  provinces  of  Earth.  Together  the 
profiles  constitute  a  three-dimensional  model  which  incorporates  considerable  structmal  detail, 
but  is  described  by  only  90  parameters  and  has  only  about  22  degrees  of  freedom.  This  is 
achieved  by  irregularly  sampling  a  detailed  regionafization  of  the  globe,  by  detailed  mapping  of 
subducted  lithosphere  in  the  mantle  as  defined  by  seismicity,  and  by  combining  these  structures 
in  an  irregular  grid  in  which  book  keeping  is  efficiently  handled.  The  resulting  RUM  model 
includes  subducting  slabs  as  sharp  fast  features  in  the  upper  mantle.  Old  continents  are  fast, 
young  oceans  are  slow.  Models  have  been  derived  for  both  compressional  and  shear  velocity.  The 
RUM  model  is  designed  to  represent  as  much  of  upper-mantle  heterogeneity  as  seen  by  body 
wave  travel  times  as  possible  with  a  simple  model.  It  can  be  useful  as  a  reference  model  for 
individual  tectonic  regions.  Travel  times  are  efficiently  generated  for  the  RUM  model. 
Mislocations  of  explosions  of  known  location  are  significantly  reduced  when  corrections  for  the 
RUM  model  are  applied  to  travel-time  residuals  for  a  spherically  symmetrical  Earth  model. 
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1.  Introduction 

Global  models  of  seismic  velocity  heterogeneity 
have  emerged  in  the  past  decade  which  show  signif¬ 
icant  deviations  from  spherical  symmetry  in  Earth’s 
structure.  Most  of  the  models  for  the  upper  mantle 
are  primarily  derived  from  surface  waves  and  yield  in¬ 
formation  about  shear- velocity  structure  [e.g.  Wood- 
house  and  Dziewonski,  1984,  Zhang  and  Tanimoto, 
1993].  As  the  quality  of  data  and  their  coverage  in¬ 
creases  more  detail  is  included  in  the  models.  Global 
models  currently  claim  to  resolve  structures  as  small 
as  1000  km  in  lateral  extent  [Zhang  and  Lay^  1996], 
although  some  authors  give  a  more  conservative  reso¬ 
lution  estimate  [  Tramperi  and  Woodhouse,  1995].  Re¬ 
gional  models,  on  the  other  hand,  resolve  structures 
to  considerably  smaller  scales,  i.e.  hundreds  of  kilo¬ 
metres  [e.g.  Zielhuis  and  Nolet,  1994,  Dehayle  and 
Leveque,  1997].  The  strength  of  heterogeneity  in  the 
models  increases  as  more  detail  is  included.  The 
global  model  RG5.5  of  Zhang  and  Tanimoto  [1993] 
has  about  ±3%  variation  at  110  km  depth  while  the 
regional  model  of  Zielhuis  and  Nolet  [1994]  has  about 
±6%  variation  at  the  same  depth. 

The  strength  of  heterogeneity  in  global  models  falls 
off  rapidly  below  a  few  hundred  km  depth.  Typically 
the  heterogeneity  is  d=3%  around  100  km  depth,  but 
less  than  ±1%  in  the  transition  zone.  Using  travel 
times,  which  potentially  offer  higher  resolution,  Gud- 
mundsson  et  al.  [1990]  found  that  the  level  of  hetero¬ 
geneity  at  large  and  intermediate  scales  (scale  >400 
km)  drops  by  about  a  factor  of  two  at  a  depth  of  300 
km.  They  also  concluded  that  smaller-scale  hetero¬ 
geneity  is  concentrated  even  more  strongly  near  the 
Earth’s  surface. 

One  important  source  of  information  about  the 
structure  of  the  upper  mantle  are  regional  refraction 
profiles  (of  length  >  2500  km)  and  body- wave  travel 
times  and  waveforms  from  earthquakes  which  can  be 
aligned  along  a  profile  of  comparable  length  and  in¬ 
terpreted  as  a  refraction  profile  [e.g.  Helmberger  and 
Wiggins,  1971,  Grand  and  Helmberger,  1984,  Walck, 
1984,  Mechie  et  al,  1993,  Kennett  et  al,  1994].  A 
comprehensive  summary  of  the  one- dimensional  mod¬ 
els  which  result  from  this  approach  can  be  found  in 
Nolet  et  al  [1994].  After  comparing  a  large  num¬ 
ber  of  models  from  varied  tectonic  environments  they 
conclude  that  there  is  clear  evidence  for  lateral  hetero¬ 
geneity  which  correlates  with  surface  tectonics.  The 
model  differences  persist  to  300-400  km  depth.  They 
also  point  out  that  local  studies  generally  produce 


larger  velocity  anomalies  than  are  found  in  global 
surface- wave  models.  This  suggests  that  global  mod¬ 
els  provide  only  a  partially  complete  picture  of  the  up¬ 
per  mantle,  with  intermediate-scale  features  filtered 
out  and  larger-scale  features  damped. 

Regional  tomography  based  on  travel  times  has  re¬ 
vealed  important  aspects  of  subducted  slab  structures 
in  the  mantle  [e.g.  Spakman  et  al,  1988,  Zhou  and 
Clayton,  1990,  van  der  Hilst  et  al,  1991,  Zhao  et  al, 
1992].  Slabs  appear  as  fast  anomalies  which  display 
a  variation  of  morphology,  particularly  at  depth.  In 
some  instances  evidence  for  necking  and  detachment 
is  seen,  some  slabs  appear  to  be  driven  through  the 
phase  transition  at  660  km  depth  while  others  pile 
up  in  the  transition  zone.  While  it  is  clear  that  slab 
morphology  is  complex  some  of  the  conclusions  from 
regional  delay-time  tomography  would  be  strength¬ 
ened  by  rigorous  hypothesis  testing.  The  most  reli¬ 
able  slab  images  are  based  on  compressional  waves. 
The  strength  of  the  anomalies  for  compressional  ve¬ 
locity  is  as  high  as  ±5%,  fast  velocities  defining  the 
slab  while  the  slower  velocities  are  generally  found  in 
the  shallow  mantle  wedge. 

The  picture  that  has  emerged  of  the  upper  man¬ 
tle  has  shear  velocity  vary  by  ±3%  on  large  scales  (> 
1000  km)  and  another  ±3%  on  intermediate  scales 
(400  -  1000  km).  Significant  heterogeneity  persists  at 
scales  smaller  yet.  Variations  of  compressional  ve¬ 
locity  are  less  well  known,  but  axe  probably  com¬ 
parable  while  somewhat  less  than  the  variations  of 
shear  velocity.  Subducting  lithosphere  possesses  an 
anomaly  of  about  5%.  Clearly,  such  strong  hetero¬ 
geneity  will  significantly  affect  the  travel  times  of 
body- wave  phases,  including  those  most  commonly 
used  in  earthquake  location.  Accounting  for  the  ef¬ 
fect  of  lateral  heterogeneity  can  therefore  significantly 
improve  the  quality  of  event  location.  Fujita  et  al. 
[1981]  estimate  that  the  structure  of  a  slab  may  affect 
the  location  of  earthquakes  within  it  by  as  much  as 
40  km.  For  teleseismic  P  waves  that  translates  into  a 
time  effect  of  about  4  seconds. 

It  is  evident  from  even  the  first  generation  of 
global  models  and  from  a  comparison  of  early  one¬ 
dimensional  models  from  different  tectonic  regions 
that  the  large-scale  features  of  shallow  upper-mantle 
structure  correlate  well  with  surface  tectonics  [Ro- 
manowicz,  1991,  Nolet  et  al,  1994].  Old  continental 
cores  have  deep  fast  roots  extending  into  the  upper 
mantle  and  spreading  centers  are  underlain  by  slow 
mantle  rocks,  although  the  depth  extent  is  subject  to 
debate.  Nataf  and  Ricard  [1996]  attach  much  signifi- 
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cance  to  the  correlation  of  tomographic  models  for  the 
upper  mantle  with  surface  tectonics  when  they  sug¬ 
gest  a  velocity  model  (3SMAC)  for  the  upper  mantle 
based  on  a  tectonic,  crustal  thickness  and  sedimentary 
thickness  regionalization,  geophysical  modeling  and 
constraints  from  other  fields  of  earth  sciences  than 
seismology.  Their  approach  is  justified  by  Ricard  et 
al.  [1996]  who  show  that  the  3SMAC  model  accounts 
well  for  normal-mode  frequencies  and  phase- velocity 
observations  world  wide. 

Like  Nataf  and  Ricard  [1996]  we  regard  the  strong 
correlation  of  surface  tectonics  and  velocity  hetero¬ 
geneity  at  shallow  depths  in  the  mantle,  where  the 
heterogeneity  is  strongest,  as  a  potential  constraint  on 
global  tomography.  In  this  paper  we  apply  a  regional¬ 
ization  as  a  constraint  to  upper-mantle  structure  de¬ 
rived  from  travel  times.  We  seek  a  simple  model  (few 
model  parameters)  to  account  for  as  much  of  upper- 
mantle  heterogeneity  as  possible.  Our  objective  is  to 
build  a  regionalized  reference  model  for  travel- time 
corrections. 

2.  Strategy 

Accounting  for  the  effects  of  velocity  heterogeneity 
on  travel  times  and  thus  event  location  must  include 
the  main  features  of  the  velocity  structure  and  re¬ 
main  a  simple  procedure  in  order  to  be  efficient.  We 
attempt  to  meet  these  prerequisites  with  the  following 
strategy: 

We  use  the  high  level  of  correlation  of  global  veloc¬ 
ity  models  with  surface  tectonics  to  define  the  main 
features  of  Earth  structure.  A  detailed  regionalization 
defines  irregularly  shaped  bodies  on  the  Earth’s  sur¬ 
face  which  in  turn  define  discrete  model  parameters. 
With  irregular  sampling  of  the  regionalization  we  can 
represent  irregular  bodies  with  a  large  rajige  of  length 
scales.  The  sampling  is  dense  where  needed  (i.e.  at 
continental  margins,  island  arcs  and  within  tectonic 
continents),  and  sparse  elsewhere  (e.g.  intraregional, 
oceanic). 

Delaunay  and  Voronoi  tesselation  are  recipes  for 
connecting  nodes  in  an  irregular  grid.  In  both  cases 
the  tesselation  is  uniquely  defined  and  automatically 
generated  given  an  arbitrary  set  of  nodes.  The  Delau¬ 
nay  tesselation  generates  a  space-filling  set  of  tetra- 
hedra  with  nodes  at  their  vertices.  The  Voronoi  tes¬ 
selation  defines  a  polyhedron  around  each  node  with 
faces  on  the  median  bisectors  between  the  node  in 
question  and  its  natural  neighbours.  [Okabe  et  al., 
1992].  Delaunay  tetrahedra  are  suitable  to  define  the 


volume  of  an  object  which  surface  has  been  sampled. 
Voronoi  polyhedra  are  suitable  to  define  the  volume 
associated  with  a  discrete  sample  of  a  varying  field. 
We  use  a  combination  of  a  Delaunay  tesselation  and 
a  Voronoi  tesselation  to  define  discrete  regions  be¬ 
longing  to  the  various  tectonic  types  [Sambridge  and 
Gudmundsson,  1997].  This  allows  us  to  significantly 
reduce  the  number  of  parameters  needed  to  represent 
the  complex  regions  compared  to  a  regular  grid  for 
example.  Efficient  algorithms  are  available  for  navi¬ 
gation  and  book  keeping  in  irregular  grids  composed 
of  Delaunay  tetrahedra  or  three-dimensional  Voronoi 
cells  [Sambridge  et  al,  1995]. 

We  use  a  select  data  set  of  travel  times  from  well 
located  events  and  nuclear  explosions  in  order  to  alle¬ 
viate  the  mapping  of  mislocation  into  structure.  We 
invert  these  data  for  an  irregularly  parameterised,  re- 
gionalised  upper-mantle  (RUM)  model.  Our  objec¬ 
tive  is  not  to  find  the  best  possible  model,  but  to  find 
the  simplest  model  that  still  contains  enough  infor¬ 
mation  about  structure  to  be  useful  for  generating 
corrections  to  travel  times  used  in  event  location. 

Teleseismic  travel  times  can  then  be  corrected  for 
this  model  by  a  relatively  simple  and  efficient  ray¬ 
tracing  procedure. 

3.  Method 

3.1.  Regionalization 

Figure  1  shows  the  regionalisation.  It  is  based 
on  maps  published  by  Sclater  et  al.  [1980]  and  Jor¬ 
dan’s  [1981]  GRTl  regionalization.  Eight  tectonic  re¬ 
gions  are  defined.  The  ocean  is  divided  into  three  age 
provinces,  the  same  as  the  GRTl  regionalization,  ex¬ 
cept  at  a  resolution  of  two  degrees  compared  to  five 
degrees  in  GRTl.  The  continents  are  divided  into 
5  regions  according  to  age,  four  regions  of  low  vol¬ 
canic  and  earthquake  activity,  in  addition  to  a  region 
termed  tectonic  continent  where  seismic  and  volcanic 
activity  is  high.  We  have  termed  the  regions  Young 
Continent,  Intermediate- age  Continent,  Old  Conti¬ 
nent  and  Ancient  Continent  as  their  separation  in 
time  does  not  strictly  correspond  to  boundaries  be¬ 
tween  geological  eras.  The  age  boundaries  are  listed 
in  Figure  1  and  were  selected  by  comparing  the  age 
zonation  of  the  continents  with  tomographic  images 
of  the  upper  mantle  [e.g.  Woodhouse  and  Dziewon- 
ski,  1984].  The  regionalization  is  very  similar  to  the 
regionalization  used  by  Nataf  and  Ricard  1996  to  gen¬ 
erate  the  3SMAC  model.  We  have  perhaps  unconven¬ 
tionally  termed  the  Iceland- Faeroe- Greenland  Ridge 
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continental  because  of  its  highly  anomalous  crustal 
thickness  for  oceanic  crust  (25-35  km  according  to 
Bott  and  Gunnarsson  [1980].  We  have  also  defined 
the  ridges  of  the  Tasman  Sea  as  continental.  A  future 
improvement  of  the  regionalization  would  include  all 
oceanic  ridges  and  plateaus  as  a  separate  oceanic  re- 
gion(s).  This  regionalization  has  obvious  limitations. 
Parts  of  the  age  mapping  involve  significant  interpo¬ 
lation  of  geo  chronological  data.  Details  are  left  out 
where  limited  by  the  resolution.  For  the  pmposes  of 
global  seismology  much  more  detail  is  not  justified 
because  teleseismic  delay  times  or  surface- wave  dis¬ 
persion  will  not  be  able  to  resolve  it. 

The  regionalization  is  initially  represented  by  16200 
parameters.  Roughly  a  quarter  of  these  are  then  se¬ 
lected  to  parameterise  the  boundaries  between  the  re¬ 
gions  (4100  parameters)  (see  Figure  lb).  The  bound¬ 
aries  are  all  parameterised  with  a  two-degree  resolu¬ 
tion  except  the  boundaries  between  oceanic  regions. 
Thus,  the  density  of  parameterisation  is  high  where 
we  expect  potentially  sharp  contrasts,  i.e.  at  conti¬ 
nental  margins  and  boundaries  between  tectonic  and 
atectonic  regions,  and  low  where  we  don’t,  i.e.  within 
the  oceans.  Figure  lb  shows  the  boundaries  of  the 
Voronoi  cells  around  the  selected  samples  of  the  re¬ 
gionalization.  The  figure  contains  about  4100  Voronoi 
cells.  Figure  Ic  shows  how  the  original  regionalization 
is  reproduced  by  the  4100  selected  cells.  The  full  com¬ 
plexity  of  Figure  la  is  reproduced  in  Figure  Ic  with 
only  about  a  quarter  of  the  number  of  parameters  in 
the  original.  Each  Voronoi  cell  is  taken  to  encompass 
a  region  of  uniform  velocity  at  each  depth.  The  re¬ 
gionalization  is  extended  from  the  surface  to  a  depth 
of  660  km  to  span  the  whole  of  the  upper  mantle. 

This  regionalization  is  not  taken  to  be  a  com¬ 
plete  description  of  velocity  heterogeneity.  We  ignore 
lower-mantle  heterogeneity  because  extrapolation  of 
surface  tectonics  into  the  lower  mantle  would  not  be 
warranted  and  because  the  strenght  of  heterogeneity 
decays  rapidly  with  depth  in  the  upper  mantle.  Shal¬ 
low  small-scale  structure  internal  to  tectonic  regions 
is  not  accounted  for  by  this  regionalisation.  In  island 
arc  regions  such  structures  are  modelled  as  subduct¬ 
ing  slabs.  In  continental  regions  such  structures  are 
accounted  for  by  station  corrections. 

3.2.  Subducting  slabs 

Subducting  slabs  are  significant  structures  in  the 
upper  mantle  which  are  clearly  oblique  and  cannot 
be  included  in  a  two-dimensional  regionalization  as 
above.  Most  of  Earth’s  seismicity  occurs  in  or  around 


subducting  slabs.  Slabs  do  affect  travel  times  signifi¬ 
cantly  and  their  presence  can,  if  not  accounted  for, 
cause  a  systematic  mislocation  of  island-arc  earth¬ 
quakes  by  as  much  as  40  km  [Fujita  et  al,  1981].  We 
have  included  oblique  slab  structures  in  our  param¬ 
eterization  of  the  upper  mantle  as  the  ninth  region. 
First  we  contoured  slab-related  seismicity  world  wide. 
A  total  of  24  separate  slab  bodies  were  defined  (see 
Figure  2a).  Examples  of  slab  contoures  are  presented 
in  Figure  2b  for  the  slabs  of  the  northwest  Pacific. 
The  contours  are  drawn  at  fifty  km  depth  intervals 
and  such  that  they  lie  near  the  top  of  the  seismo- 
genic  region.  We  used  the  relocated  catalog  of  Eng- 
dahl  et  al  [1997]  as  well  as  a  filtered  version  of  the 
ISC  catalog  (based  on  location  quality)  to  define  the 
seismicity.  Neither  catalog  is  complete  and  both  con¬ 
tain  significant  gaps  which  are  difficult  to  interpret. 
These  gaps  were  not  interpolated  if  they  exceed  a  few 
hundred  km  in  size.  The  slab  contours  obviously  do 
not  define  a  slab  where  there  is  no  seismicity.  It  is  dif¬ 
ficult  to  invisage  the  thermal  anomaly  of  subducting 
lithosphere  to  terminate  as  abruptly  as  the  seismic¬ 
ity  in  many  places.  Thus  our  model  of  slabs  is  likely 
to  be  incomplete.  However,  slabs  are  defined  where 
we  most  need  them,  i.e.  where  the  earthquakes  occur 
and  for  which  travel- time  corrections  are  needed. 

The  slab  contours  define  a  surface  which  we  take 
to  define  the  top  of  the  slab  (or  what  used  to  be  the 
top  of  the  oceanic  lithosphere  before  subduction).  We 
then  define  a  complimentary  surface  assuming  that 
the  slab  is  200  km  thick.  This  is  a  simplification  of 
slab  geometry,  particularly  at  depth.  No  necking  or 
slab  detachment  is  allowed  for,  while  it  is  possible 
that  that  type  of  morphology  is  associated  with  gaps 
in  seismicity  [Widiyantoro  and  van  der  Hilst,  1996]. 
Slabs  are  not  allowed  to  pile  up  in  the  transition  zone 
as  appears  to  be  the  case  for  e.g.  the  Izu-Bonin  slab 
[van  der  Hilst  et  a/.,  1993]. 

Projecting  the  original  surface  200  km  along  its  lo¬ 
cal  normal  occasionally  results  in  the  folding  of  the 
original  surface.  Such  folds  were  simply  edited  out 
interactively  with  a  graphics  tool.  We  sample  the 
contours  using  nodes  placed  at  a  lateral  spacing  of 
about  10  km.  These  nodes  are  not  internal  to  the 
slabs  as  we  define  them,  but  on  their  surface.  In 
this  case  const  ant- velocity  Voronoi  cells  as  were  used 
in  the  surface  parameterization  are  not  suitable  to 
describe  their  complex  geometry,  because  these  cells 
surround  each  node  and  would  necessarily  extend  the 
slab  and  change  its  volume.  Instead  we  place  the 
nodes  of  the  slabs  in  a  global  coordinate  system  and 
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construct  the  three-dimensional  Delaunay  tesselation 
for  the  approximately  23,000  slab  samples  for  all  24 
individual  slab  structures.  Those  Delaunay  tretra- 
hedra  which  join  the  two  surfaces  of  the  same  slab 
structure  then  define  the  volume  of  slabs  in  the  up¬ 
per  mantle.  The  Delaunay  tetrahedral  representation 
of  the  slabs  of  the  northwest  Pacific  is  shown  in  per¬ 
spective  in  Figure  2b  as  an  example.  These  are  some 
of  the  larger  and  deeper  slabs  which  display  signifi¬ 
cant  complexity  in  shape. 

3.3.  Summary  of  parameterization 

In  summary  the  parameterization  of  the  upper 
mantle  consists  of  two  distinct  elements:  Delaunay 
tetrahedra  to  define  slabs  and  a  two-dimensional  grid 
of  Voronoi  cells  extended  uniformly  to  a  depth  of  660 
km  to  define  tectonic  regions.  The  slabs  are  overlaid 
on  the  regionalization.  To  determine  which  type  of 
environment  a  given  point  (along  a  ray)  is  in  we  first 
detect  if  it  is  within  a  slab.  If  it  is  not,  its  position 
in  the  regionalization  is  determined.  Thus  we  repre¬ 
sent  complex  structures  efficiently  and  economic  algo¬ 
rithms  are  available  to  locate  any  arbitrary  point  in 
either  irregular  grid  [see  Sambridge  and  Gudmunds- 
son,  1997].  The  parameterisation  includes  fine  spatial 
features  in  subduction  zones  on  a  scale  of  10  km  and 
also  tectonic  age  provinces  a  few  thousand  kilometers 
across.  Its  dynamic  range  is  thus  two  to  three  orders 
of  magnitude. 

3.4.  Data  selection 

We  use  a  select  data  set  of  travel  times  from  well  lo¬ 
cated  events  and  nuclear  explosions  in  order  to  allevi¬ 
ate  the  mapping  of  mislocation  into  structure.  This  is 
of  particular  importance  in  building  a  reference  Earth 
model  because  observed  biases  in  earthquake  location 
can  be  significant  and  can  inject  artificial  structural 
signal  into  travel  times  on  the  order  of  seconds.  The 
events  we  use  are  the  105  events  of  the  test-data  set 
used  by  Kennett  and  Engdahl  [1991]]  supplemented 
with  92  large  earthquakes  to  get  a  more  even  sam¬ 
pling  of  the  globe  and  the  9  tectonic  regions.  The 
criterion  for  the  selection  of  these  events  was  that 
they  be  recorded  by  a  large  number  of  globally  dis¬ 
tributed  stations.  Their  locations  are  taken  from  nu¬ 
clear  agencies,  local  and  regional  networks,  and  if  nei¬ 
ther  of  the  above  is  available  from  the  relocated  cat¬ 
alog  of  Engdahl  et  al  [1997].  They  reassociated  ISC 
data,  included  pP  and  pwP  phases  in  the  location, 
used  IASP91  [Kennett  and  Engdahl,  1991]  and  sub¬ 
jected  locations  to  stringent  fitness  criteria.  We  use 


only  teleseismic  P  and  PKIKP  arrivals  (first  arrivals) 
for  the  compressional- velocity  model  and  teleseismic 
S  arrivals  for  the  shear-velocity  structure.  We  use 
only  picks  reported  to  the  ISC  as  impulsive  to  insure 
a  high  degree  of  robustness  of  individual  picks. 

All  the  data  have  been  corrected  for  the  Earth’s 
elliptic! ty  [Kennett  and  Gudmundsson,  1996],  for  sta¬ 
tion  elevation,  and  for  surface  topography /bathymetry 
in  the  case  of  phases  reflecting  off  the  Earth’s  surface. 
We  have  prepared  data  from  other  phases  such  as  pP, 
PP,  PcP,  PKP,  sS,  SS,  ScS,  and  SKS  from  the  197 
events  in  this  select  data  set. 

Figure  3  shows  a  map  of  the  197  events  used 
together  with  the  regionalization.  The  events  are 
grouped  into  three  categories.  First  there  are  21 
nuclear  and  chemical  explosions  (displayed  as  white 
stars  in  the  figure).  Their  locations  will  be  taken  to 
be  fixed.  Then  there  are  39  earthquakes  which  are 
well  located  with  local  observations  in  a  local  veloc¬ 
ity  model.  They  are  displayed  as  white  dots.  Their 
locations  are  also  assumed  fixed.  The  remaining  137 
events  (shown  as  black  dots)  are  located  using  tele¬ 
seismic  observations  and  are  allowed  to  vary  in  their 
origin  time  and  depth. 

We  reference  travel-time  residuals  to  the  AK135 
model  [Kennett  et  al,  1995]  and  reassociate  all  im¬ 
pulsive  picks  reported  to  the  ISC  for  the  above  events. 
We  find  a  surprising  number  of  events  with  a  signifi¬ 
cant  mean  residual.  This  is  surprising  because  AK135 
is  almost  identical  in  the  mantle  to  IASP91  which  En¬ 
gdahl  et  al.  [1997]  used  for  their  locations.  The  ex¬ 
planation  must  lie  in  a  difference  in  data  selection. 
We  for  example  include  PKIKP  and  not  pP.  We  do 
not  associate  pwP.  The  large  number  of  events  with  a 
significant  apparent  origin-time  shift  led  us  to  allow 
origin  time  to  be  a  free  parameter  in  our  inversion  for 
regionalized  structure.  Allowing  origin  time  to  vary 
independently  for  P  and  S  is  equivalent  to  allowing  for 
a  relocation  in  source  depth  and  origin  time.  We  find 
14291  P  and  PKIKP  picks  and  2608  S  picks  which 
satisfy  our  selection  criteria, 

4.  Inversion  -  data  reduction 

We  solve  a  linearized  inverse  problem.  A  higher 
order  approach  is  not  warranted  because  ray  bend¬ 
ing  critically  depends  on  velocity  gradients  which  are 
artificially  represented  by  our  parameterization.  We 
sample  depth  ten  times  at  66  km  spacing  and  seek 
9  independent  velocity  profiles,  one  for  each  of  the 
eight  tectonic  regions  as  well  as  slabs.  We  use  a 
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standard  damped-least-squares  inversion  procedure. 
In  regions  with  limited  depth  distribution  of  sources 
teieseismic  rays  give  limited  depth  resolution  in  the 
upper  mantle.  We  cannot  well  resolve  ten  indepen¬ 
dent  parameters  except  in  subducting  slabs,  where 
there  is  a  distribution  of  events  in  depth.  Since  the 
level  of  heterogeneity  in  the  upper  mantle  is  strongly 
concentrated  near  the  surface  [e.g.  Gudmundsson  et 
al,  1990],  in  the  top  300-400  km,  we  have  projected 
the  depth  profiles  on  to  a  set  of  four  Gaussian  ba¬ 
sis  functions  which  are  not  allowed  to  put  structure 
in  the  lower  half  of  the  upper  mantle.  In  effect  we 
are  implementing  depth-dependent  damping.  We  ar¬ 
gue  that  this  is  in  line  with  our  general  approach  of 
applying  stringent  geometrical  constraints  to  the  in¬ 
version  based  on  a  priori  information.  The  basis  func¬ 
tions  are  not  orthogonal  in  themselves  but  decompose 
into  a  fourfold  orthogonal  basis.  The  model  for  sub¬ 
ducted  lithosphere  is  allowed  to  vary  throughout  the 
upper  mantle.  The  results  we  present  are  moderately 
damped.  The  trade  off  employed  to  select  the  damp¬ 
ing  was  that  of  variance  reduction  versus  damping 
coefficient  (which  varies  monotonically  with  resolu¬ 
tion).  The  damping  was  selected  where  it  started  to 
significantly  affect  the  variance  reduction  (lowered  by 
10  —  20%).  The  relocation  component  of  the  inverse 
problem  is  scaled  to  allow  for  virtually  full  removal  of 
average  residuals  for  each  event. 

For  the  P-wave  data  origin-time  shifts  of  137  events 
out  of  192  account  for  about  55%  variance  reduc¬ 
tion  (see  Table  1).  An  undamped  inversion  allowing 
for  10  independent  depth  parameters  in  each  region 
achieves  a  19%  variance  reduction.  Limiting  the  de¬ 
grees  of  freedom  to  four  per  region  via  the  Gaussian 
basis  reduces  the  variance  reduction  to  16%.  We  have 
damped  the  solution  to  yield  a  variance  reduction  of 
14%.  For  the  S-wave  data  the  pattern  is  similar,  with 
the  damped  model  achieving  a  10%  variance  reduc¬ 
tion.  For  the  S-waves  the  variance  is  considerably 
higher  than  for  the  P-waves  and  the  relative  variance 
reduction  achieved  by  origin  time  shift  is  much  less 
than  for  the  P-wave  data. 

We  find  that  none  of  the  minor  data  sets  consisting 
of  the  compressional  pP,  PP,  PcP,  and  PKP  phases 
or  the  shear  sS,  SS,  ScS,  and  SKS  phases  improve 
resolution  significantly  nor  are  they  consistent  with 
the  major  data  sets  of  P,  PKIKP,  or  S  phases.  The 
RUM  model  does  not  significantly  reduce  the  variance 
of  the  minor  data  sets. 

As  a  final  step  in  the  data  reduction  we  estimate 
station  corrections  to  be  appHed  after  correction  for 


the  RUM  model.  We  evaluate  a  correction  only  for 
stations  with  more  than  4  picks  reported.  The  station 
correction  is  simply  the  average  residual  for  each  sta¬ 
tion  after  correction  for  the  RUM  model.  This  step 
achieves  a  further  variance  reduction  of  16%  and  11%) 
for  the  P  waves  and  S  waves,  respectively.  Table  1 
summarizes  the  reduction  of  the  mean  and  variance 
of  the  P-  and  S-wave  data  by  the  individual  steps  of 
the  data  reduction.  Figure  4  shows  a  histogram  of 
time  residuals  before  and  after  data  reduction.  We 
note  that  the  quoted  variance  reduction  applies  to  in¬ 
dividual  picks,  not  summary  rays.  We  also  note  that 
the  remaining  variance  in  the  data  is  comparable  to, 
but  slightly  higher  than  the  estimates  of  noise  vari¬ 
ance  in  the  ISC  data  for  teieseismic  P  and  S  waves  by 
Gudmundsson  et  al.  [1990]  and  Davies  et  al  [1992], 

5.  Results 

The  results  of  the  inversion  are  shown  in  Figure  5 
and  Table  2.  In  Figure  5  each  function  of  depth  rep¬ 
resents  a  velocity  perturbation  for  a  given  tectonic  re¬ 
gion.  Figure  5a  shows  the  results  for  P  waves,  Figure 
5b  for  S  waves.  The  results  are  plotted  as  a  relative 
velocity  perturbation  (in  percent)  with  one  standard 
deviation  error  bars  assuming  that  the  error  variance 
is  as  given  by  Gudmundsson  et  al.  [1990]  for  the  tele- 
seismic  P-wave  data,  i.e.  Is^,  and  by  Davies  et  al. 
[1992]  for  teieseismic  S-wave  data,  i.e  4s^,  Table  2 
also  list  the  results  as  relative  velocity  perturbations. 

The  results  for  compressional  velocity  are  encour¬ 
aging.  Slabs  appear  as  a  continuous  fast  anomaly 
throughout  the  upper  mantle  except  at  the  top  where 
we  have  tectonic  crust  rather  than  a  true  slab  and 
where  the  result  is  probably  affected  by  the  slow  man¬ 
tle  wedge.  The  amplitude  of  the  slab  structure  is 
2-3%  velocity  perturbation.  This  is  assumed  to  be 
uniform  over  a  thickness  of  200  km.  The  subduct¬ 
ing  lithosphere  is  probably  not  that  thick  at  shallow 
depths,  but  the  associated  thermal  anomaly  diffuses 
with  time  or  depth.  Furthermore,  it  would  peak  above 
its  average  level.  The  2-3%  level  of  perturbation  in 
the  slabs  would  thus  be  consistent  with  a  somewhat 
higher  peak  perturbation  in  the  core  of  the  slabs, 
possibly  as  high  as  5%.  This  is  a  reasonable  level 
of  perturbation  compared  to  predictions  on  the  ba¬ 
sis  of  thermal  modeling.  Other  features  of  the  results 
are  also  encouraging.  The  oceans  are  slow  with  old 
oceans  developing  a  fast  lithosphere.  Tectonic  conti¬ 
nent  and  young  continent  are  also  slow  and  essentially 
indestinguishable.  This  is  of  interest  because,  e.g. 
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Nataf  and  Ricard  [1996]  do  not  make  a  distinction  be¬ 
tween  those  regions  in  their  regionalization.  That  de¬ 
cision  appears  to  have  been  justified.  Currently  active 
continental  margins  and  island  arcs  are  not  clearly 
distinguishable  from  other  Mesozoic  continental  re¬ 
gions.  Intermediate-age  continents  have  developed  a 
fast  lithosphere  while  their  asthenosphere  looks  much 
like  that  of  young  continental  regions.  There  is  no  re¬ 
solvable  structure  underneath  old  continents  suggest¬ 
ing  that  AK135  (and  IASP91)  are  representative  of 
continents  in  that  age  range.  This  is  interesting  and 
not  surprising  because  the  AK135  and  IASP91  models 
were  not  designed  to  represent  average  Earth  struc¬ 
ture  but  rather  as  effective  reference  models  for  travel 
time  prediction  deliberately  biased  according  to  the 
distribution  of  stations.  Ancient  continents  are  fast, 
but  not  to  a  great  depth,  only  about  200  km.  Note, 
however,  that  depth  resolution  is  limited.  It  should 
be  noted  that  we  have  not  corrected  for  crustal  struc¬ 
ture  so  the  model  should  be  interpreted  to  include 
crustal  signature.  The  number  of  degrees  of  freedom 
in  the  compressional  model  is  22  (trace  of  resolution 
matrix). 

The  RUM  model  is  best  compared  to  individual 
one-dimensional  models  by  comparing  the  integral 
travel-time  anomaly  at  vertical  incidence.  The  RUM 
model  gives  a  range  from  -0.25  to  1.16  seconds  relative 
to  AK135.  The  negative  (fast)  extreme  corresponds 
to  ancient  continent  while  the  positive  (slow)  extreme 
corresponds  to  young  ocean.  Shield  models  such  as 
S25  of  Lefevre  and  Helmberger  [1989]  for  the  Lauren- 
tian  shield,  quartzn  of  Mechie  et  al.  [1993]  for  the 
Fennorussian  shield,  and  fish  of  Kennett  et  al.  [1994] 
for  the  Australian  shield,  have  an  average  travel-time 
anomaly  of  about  -1  second.  The  tectonic  models  t7 
and  t9  of  Burdick  and  Helmberger  [1978]  and  Burdick 
[1981]  average  at  about  0.6  seconds.  Tectonic  conti¬ 
nent  in  RUM  has  an  anomaly  of  about  0.5  seconds. 
The  oceanic  models  gca  and  cjf  of  Walck  [1984,  1985] 
average  at  0.2  seconds.  The  range  is  the  same  for 
the  elements  of  RUM  and  the  one-dimensional  mod¬ 
els,  i.e.  about  1.5  seconds.  The  time  anomaly  for 
tectonic  continent  is  similai'  in  RUM  as  in  tectonic 
one-dimensional  models,  but  for  shield  models  and 
oceanic  models  there  is  a  shift  of  about  1  second.  We 
suggest  a  few  potential  explanations  of  this  discrep¬ 
ancy.  Stations  and  events  are  co-located  in  tectonic 
regions.  Therefore,  one-dimensional  refraction  type 
studies  enjoy  good  constraints  on  shallow  structure. 
This  is  not  the  case  in  shield  regions  or  oceanic  re¬ 
gions.  The  stations  are  removed  from  the  seismicity. 


Therefore,  the  shallow  constraint  is  poor  and/or  the 
paths  tectonically  mixed  at  short  distance.  From  Fig¬ 
ure  3  it  is  evident  that  we  do  not  have  any  events 
in  young  oceanic  regions  that  are  constrained  by  lo¬ 
cal  observations.  Neither  do  we  have  such  events  in 
shield  regions.  The  RUM  model  may  therefore  suffer 
from  a  lack  of  absolute  reference  in  those  regions.  On 
the  other  hand  stations  in  shield  regions  record  many 
of  the  locally  constrained  events  which  adds  absolute 
constraint  to  the  shields.  The  RUM  model  represents 
an  average  property  for  all  regions  of  the  same  tec¬ 
tonic  type.  Station  corrections  (discussed  below)  in¬ 
dicate  that  considerable  variability  exists  among  re¬ 
gions  of  the  same  type.  It  is  possible  that  published 
one-dimensional,  regional  models  are  not  representa¬ 
tive  of  their  regions. 

The  results  for  the  shear  velocity  are  less  well  con¬ 
strained  and  more  poorly  resolved.  This  is  because 
of  fewer  and  poorer  data.  In  order  to  render  the  re¬ 
sult  interpretable  in  the  sense  of  suppressing  noise 
sufficiently  so  that  model  exceeds  uncertainty  in  am¬ 
plitude,  we  have  had  to  damp  the  S-wave  model  more 
heavily  than  the  P-wave  model.  The  number  of  de¬ 
grees  of  freedom  in  the  shear- velocity  model  is  only  1 1 . 
Still  subducted  lithosphere  extends  as  a  fast  anomaly 
of  reasonable  amplitude  throughout  the  upper  man¬ 
tle,  old  continental  cores  are  fast  and  oceans  are  slow. 
The  absence  of  a  lithosphere  in  the  oceans  may  be 
ascribed  to  lower  resolution.  The  reduced  level  of 
coherency  in  the  transition  from  young  ocean  to  an¬ 
cient  continent  can  be  understood  in  terms  of  poorer 
constraint.  The  range  of  the  predicted  vertical  travel¬ 
time  anomaly  for  RUM-S  is  -2.1  seconds  for  ancient 
continents,  which  is  comparable  to  the  sna  model  of 
Grand  and  Helmberger  [1984]  and  the  njpb  model  of 
Kennett  et  al.  [1994]  for  the  North- American  and 
Australian  shields,  respectively,  to  about  2  seconds  in 
the  oceans.  The  vertical  time  delay  for  the  tectonic 
regions  of  RUM-S  is  only  about  1  second  while  the  tna 
model  of  Grand  and  Helmberger  [1984]  has  a  7.5  sec¬ 
ond  delay.  Here  the  agreement  in  the  tectonic  region 
is  of  the  same  sense  but  poor  in  amplitude. 

We  refrain  from  interpreting  the  results  beyond 
pointing  out  that  dispite  strict  geometrical  constraints 
on  the  lateral  extent  of  anomalies  the  level  of  velocity 
perturbation  is  reasonable.  Any  comparison  between 
the  results  for  P  waves  and  S  waves  would  have  to  take 
into  account  the  difference  in  the  damping  of  the  two 
models.  The  S-wave  model  is  more  heavily  damped 
than  the  P-wave  model,  which  qualitatively  explains 
why  the  modeled  perturbation  in  shear  velocity  does 
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not  cleaxly  exceed  the  perturbation  of  compressional 
velocity  (we  would  expect  the  former  to  exceed  the 
latter  by  50  —  100%).  Surface  waves  may  be  better 
suited  to  study  the  shear- velocity  heterogeneity  of  the 
upper  mantle  than  delay  times. 

Our  confidence  in  the  results  is  much  higher  for 
compressional  velocity  than  shear  velocity.  It  is  inter¬ 
esting  to  display  the  model  in  its  full  tree-dimensionality. 
Figure  6  shows  a  pair  of  cross  sections  through  the 
RUM-P  model.  In  Figure  6a  the  section  cuts  through 
the  Tongan  slab,  the  Banda  Sea  slab  near  its  sharp 
bend,  through  some  shallow  slabs  in  the  Philippines, 
and  through  the  Andean  slab.  It  also  grazes  through 
the  Hellenic  slab.  In  Figure  6b  the  cross  section 
cuts  through  the  Sumateran  slab,  shallow  slabs  in 
the  Phillippines  and  the  slabs  of  the  northwestern  Pa¬ 
cific  along  strike.  Subducting  lithosphere  constitutes 
a  major  feature  of  the  RUM  model  at  the  extrem¬ 
ity  of  fast  perturbations.  Their  geometry  is  imposed. 
This  may  render  the  estimate  of  their  overall  veloc¬ 
ity  perturbation  more  robust  than  in  less  constrained 
delay-time  models  of  slab  regions.  The  velocity  per¬ 
turbations  in  the  upper  mantle  range  between  ±2.5%, 
which  is  similar  to  what  one  might  expect  on  the  basis 
of  estimates  of  i/  =  dlnVs/dlnVp  and  the  heterogene¬ 
ity  pattern  as  modeled  from  surface-wave  data. 

Figure  7  shows  the  station  corrections  for  P-wave 
residuals,  A  corresponding  map  for  S-wave  resid¬ 
uals  is  not  presented  because  of  its  sparseness  and 
how  poorly  constrained  we  regard  the  shear- velocity 
model.  Clearly  a  significant  signal  is  left  in  the  data 
after  inversion  for  regionalized  structure.  Note  for 
example  that  a  fast  anomaly  remains  within  the  Lau- 
rentian  craton.  On  the  other  hand  the  African  era- 
tons  are  left  with  slow  anomalies.  This  is  in  agree¬ 
ment  with  results  of  global,  surface-wave  tomogra¬ 
phy  which  has  the  African  cratons  show  up  much 
less  clearly  than  e.g.  the  Laurentian,  Fennorussian 
and  Australian  cratons.  Some  tectonic  regions  are 
left  with  systematically  fast  residuals  while  others  are 
slow.  Considerable  scatter  is  evident  within  some  re¬ 
gions.  Clearly  all  regions  of  the  same  type  are  not  the 
same.  Furthermore,  small-scale  structure  within  in¬ 
dividual  regions  is  not  modeled,  but  instead  absorbed 
in  station  corrections. 

6.  Locating  known  events 

To  test  the  success  of  the  RUM  model  as  a  ref¬ 
erence  for  teleseismic  travel  times  we  relocate  19  of 
the  21  explosions  in  the  data  set.  The  locations  of 


those  19  explosions  can  be  regarded  as  known  [Ken- 
nett  and  Engdahlf  1991].  We  compare  relocations  us¬ 
ing  teleseismic  P  residuals  which  are  referenced  to 
the  AK135  model,  to  the  RUM-P  model,  and  to  the 
RUM-P  model  as  well  as  station  corrections.  We  re¬ 
locate  only  the  epicenter,  i.e.  the  depth  is  fixed  at 
the  reported  depth.  Many  of  the  earlier  events  in 
particular  have  a  very  unevenly  distributed  sampling. 
The  locations  become  much  more  robust  by  evening 
this  distribution  or  declustering  it.  For  each  event 
we  construct  a  kernel  matrix  which  is  nx3,  where  n 
is  the  number  of  observations.  The  three  elements  of 
each  row  are  the  sensitivity  to  relocation  in  latitude, 
longitude,  and  origin  time  (the  last  is  always  -1).  We 
bin  the  first  and  second  columns  of  the  kernel  matrix 
four  times  each.  After  removing  outliers  (residuals 
more  than  5  seconds  away  from  mean  of  bin)  we  cal¬ 
culate  an  average  residual  for  each  of  the  sixteen  bins 
downweighting  emergent  picks.  This  procedure  has 
a  similar  effect  in  improving  the  robustness  of  loca¬ 
tion  to  the  weighting  scheme  employed  by  Smith  and 
Ekstrom  [1995].  We  use  only  P  residuals  in  the  tele¬ 
seismic  range.  Finally,  we  solve  a  linearized  relocation 
problem  giving  all  bins  equal  weight.  In  effect  we  are 
declustering  the  sampling  of  the  focal  sphere  as  seen 
from  the  event  to  minimize  biases  in  the  relocation. 
The  procedure  of  the  relocation  is  the  same  for  all 
events  and  for  all  levels  of  data  reduction.  The  re¬ 
sults  of  the  relocations  are  listed  in  Table  4  as  misfit 
in  origin  time  (seconds)  and  length  of  the  horizon¬ 
tal  mislocation  vector  (kilometers)  and  summarized 
at  the  bottom  as  the  root-mean-square  (rms)  misfit 
in  time  and  horizontally.  The  mislo  cations  of  indi¬ 
vidual  events  using  AK135  for  reference  are  similar 
to  the  mislo  cations  reported  by  Kennett  and  Eng- 
dahl  [1991].  In  other  words  the  level  of  robustness  of 
the  location  procedure  is  similar  although  the  proce¬ 
dures  differ.  After  correction  of  residuals  for  the  RUM 
three-dimensional  model  the  mislocations  are  signif¬ 
icantly  improved,  i.e.  by  about  20%.  Applying  the 
RUM  station  corrections  reduces  the  rms  misfit  by  a 
further  33%.  While  it  is  tempting  to  compare  the  rms 
misfit  to  that  achieved  by  Smith  and  Ekstrom  [1995] 
in  absolute  terms  that  would  have  limited  meaning 
because  the  location  procedure  is  not  the  same. 

A  significant  level  of  reduction  of  mislocation  is 
related  to  the  application  of  station  corrections.  The 
RUM  model’s  usefulness  for  travel-time  corrections 
lies  in  the  fact  that  the  corrections  are  derived  from  a 
model  and  that  it  effectively  implements  corrections 
for  near-source  structure,  event  corrections.  Where 
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a  stations  history  does  not  allow  for  a  robust  station 
correction  the  model  provides  a  reasonable  means  of 
predicting  that  station  correction  to  first  order. 

7.  Conclusions 

We  have  described  upper-mantle  structure  with  22 
degrees  of  freedom  (compressional  velocity)  and  11 
degrees  of  freedom  (shear  velocity)  and  managed  to 
explain  a  significant  portion  of  the  apparent  struc¬ 
tural  variance  in  travel-time  data.  The  data  are  indi¬ 
vidual,  not  summary-ray  data.  S-wave  residuals  give 
results  similar  to  the  P-wave  velocity  variations,  but 
less  well  constrained.  Our  parameterization  in  effect 
imposes  complex  spatial  constraints  on  this  inversion 
for  upper-mantle  structure  which  we  can  think  of  as 
a  priori  information,  i.e.  the  presumption  that  the 
upper  mantle  can  be  regionalized  as  we  have  done. 
The  success  of  the  inversion  justifies  those  constraints. 
The  remaining  station  residuals  point  to  the  fact  that 
this  regionalized  description  is  not  complete,  however, 
and  that  on  order  of  half  the  structural  signal  in  the 
data  is  explained  by  it.  The  RUM  model  is  derived 
as  a  regionalized  reference  model  for  the  upper  man¬ 
tle.  It  should  prove  useful  for  correcting  travel  times 
for  upper-mantle  heterogeneity  to  first  order.  Im¬ 
provements  in  the  location  misfit  of  known  explosions 
demonstrates  the  models  success.  The  RUM  model 
can  also  be  used  as  a  starting  point  for  tomography, 
an  initial  starting  model. 

The  RUM  model  is  tabulated  in  Table  2.  We  have 
also  prepared  a  home  page  on  the  world-wide  web 
at  http://rses.anu.edu.au/RUM/rum.html  containing 
graphical  information,  software  and  software  descrip¬ 
tion.  The  software  includes  codes  to  calculate  travel¬ 
time  corrections  for  the  RUM  model  for  a  variety  of 
phases  including  other  first  order  effects  such  as  due 
to  the  Earth’s  topography. 
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Figure  1.  Tectonic  regionalization  and  its  reproduction  using  Voronoi  cells,  a)  the 
regionalization  as  designed  on  a  2x2  degree  grid  using  16200  parameters,  b)  boun¬ 
daries  of  Voronoi  cells  around  4100  nodes  selected  to  sample  the  details  of  the  region¬ 
alization,  c)  the  regionalization  as  represented  by  4100  Voronoi  cells. 

Figure  2.  a)  Contoured  subducting  slabs  with  topographic/bathymetric  map  shown  for 
reference,  b)  Contours  of  slabs  of  the  northwest  Pacific  with  topographicA)athymetric 
map  for  reference,  c)  A  perspective  view  of  the  slabs  of  the  northwest  Pacific  from 
above  and  the  southwest.  The  edges  of  individual  Delaunay  tetrahedra  are  shown  in 
plack. 

Figure  3.  Epicenters  of  the  197  events  used  in  this  study  with  regionalization  for 
reference.  White  stars  represent  artificial  explosions,  white  dots  represent  events 
located  within  a  regional  network  of  stations,  and  black  dots  represent  events  located 
using  global  observations. 

Figure  4.  Histograms  show  distribution  of  time  residuals  before  and  after  inversion/ 
data  reduction,  a)  P-wave  residuals,  b)  S-wave  residuals. 

Figure  5.  The  RUM  velocity  model  presented  as  percentile  velocity  perturbation  from 
the  AK135  model,  a)  For  compressional  velocity,  b)  For  shear  velocity.  The  shaded 
area  represents  the  range  ±a,  where  a  is  one  standard  deviation. 

Figure  6.  Two  cross  sections  through  the  three-dimensional  RUM  model  showing  the 
range  of  velocity  perturbations  and  demonstrating  the  types  of  bodies  included  in  the 
model.  Note  that  this  global  model  contains  clear  slabs  in  the  upper  mantle. 

Figure  7.  Station  corrections  for  P-wave  residuals  after  correction  for  the  RUM 
model.  White  squares  represent  negative  corrections  (fast),  black  triangles  represent 
positive  corrections  (slow).  The  symbol  area  is  proportional  to  the  amplitude  of  the 
correction  with  1  second  corrections  given  in  legend  for  reference. 

Table  1.  Mean  and  variance  reduction  after  the  different  steps  of  inversion/data  reduc¬ 
tion.  Upper  for  P  residuals,  lower  for  S  residuals. 

Table  2.  The  RUM  model  as  a  percentile  slowness  perturbation  relative  to  AK135. 
The  labeling  of  the  regions  is  as  in  Figure  5.  Upper  for  compressional  velocity,  lower 
for  shear  velocity. 

Table  3.  Mislocation  of  explosions  of  known  location  after  correction  for  the  RUM 
model  and  the  associated  station  corrections.  Only  teleseismic  residuals  were  used. 
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